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ABSTRACT 


The  attitude  dynamics  of  a  gravity -gradient  oriented  satellite,  which  employs  soft  nickel- 
iron  rods  as  magnetic  hysteresis  dampers,  is  simulated  by  a  digital  computer  program. 

A  subroutine  was  generated  to  compute  the  flux  induced  in  the  rods  for  an  arbitrary  ve¬ 
hicle  orientation  relative  to  the  earth's  magnetic  field  vector.  This  subroutine  was 
added  to  an  existing  computer  program.  Analytical  expressions  were  fitted  to  empirical 
major  and  minor  hysteresi3  loops,  and  a  logic  scheme  was  devised  for  tracing  BH  curves 
with  an  arbitrarily  fluctuating  applied  amgnetic  field.  This  enables  simulation  of  the 
effects  due  to  the  magnetic  field  encountered  in  orbit.  Volume  I  of  this  report  describes 
the  analysis  and  the  computer  program  in  detail.  Volume  EL  presents  the  results  of 
computer  runs  made  under  various  initial  conditions. 
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SECTION  1. 

INTRODUCTION 

The  U.  S.  Naval  Weapons  Laboratory  at  Dahlgren,  Va. ,  awarded  the  Spacecraft  Department 
of  the  General  Electric  Co.  a  contract  for  performing  the  digital  simulation  of  the  attitude 
dynamics  of  a  gravity-gradient  oriented  satellite  employing  soft  nickel-iron  rods  as  mag¬ 
netic  hysteresis  dampers.  The  work,  performed  under  Contract  No.  N178-8450  as 
amended  by  Contract  Change  Notice,  consisted  of  (1)  obtaining  test  data  for  the  rod  mate¬ 
rial  from  the  Allegheny  Ludlum  Steel  Corporation;  (2)  fitting  analytical  expressions  to  these 
magnetic  curves;  (3)  generating  a  magnetic  torque  subroutine  for  incorporation  into  an 
existing  GE  digital  program  for  simulating  the  attitude  dynamics  of  the  satellite;  and  (4) 
making  computer  runs  with  various  initial  conditions  and  vehicle  parameters,  as  specified 
by  NWL. 

The  report  on  this  contract  is  submitted  in  two  volumes.  Volume  I  includes  a  summary  of  the 
work  done,  explains  the  magnetics  fundamentals  involved,  and  describes  the  approaches  taken 
to  solve  the  problems.  A  complete  description  of  the  hysteresis  torque  subroutine  is  pre¬ 
sented,  including  the  logical  rules,  the  scheme  for  implementing  the  logic,  and  the  complete 
equations.  The  unsuccessful  approaches  are  also  described  briefly. 

Volume  II  includes  a  brief  description  of  the  satellite  in  terms  of  the  characteristics  and 
parameters  which  affect  its  attitude  performance,  a  list  of  the  initial  and  other  conditions 
for  the  various  computer  runs,  and  the  results  of  these  runs. 

The  Allegheny  Ludlum  data  used  is  for  the  same  material  as  that  in  the  hysteresis  rods  in 
the  satellite.  However,  the  heat  treatment  of  the  satellite  rods  is  slightly  different  from 
that  of  the  sample  used  for  taking  the  data.  Also,  the  maximum  value  of  the  earth's  magnetic 
field  encountered  in  orbit  is  somewhat  greater  than  the  maximum  value  used  in  taking  the 
data.  The  data  were  therefore  extrapolated  to  provide  the  required  range.  Additional  data, 
taken  with  greater  field  strengths  applied  to  samples  having  appropriate  heat  treatment, 
have  been  requested  from  Allegheny  Ludlum,  but  have  not  yet  been  received.  When  they 
are  received,  these  curves  can  be  compared  with  those  used  in  the  digital  simulation. 
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SECTION  2. 

SUMMARY 

The  work  reported  herein  includes  (1)  determining  the  general  requirements  of  the  hysteresis 
torque  subroutine,  in  terms  of  the  representation  of  the  magnetic  curves  and  their  use  in 
conjunction  with  a  logic  scheme  in  a  m?nner  which  is  compatible  with  the  digital  computer; 

(2)  fitting  analytical  expressions  to  the  original  magnetic  data,  after  making  the  corrections 
for  the  demagnetization  factor  of  the  rods;  ana  (3)  devising  a  scheme  for  using  these  ex¬ 
pressions,  together  with  a  minimum  amount  of  information  about  the  past  history  of  each 
rod,  to  derive  magnetic  curves  under  the  conditions  associated  with  an  arbitrarily  varying 
applied  magnetic  field. 

The  magnetics  and  dynamics  fundamentals  which  bear  on  the  problem  are  discussed.  It  is 
shown  that  a  continuously  varying  magnetic  field  applied  to  a  hysteresis  material  causes  a 
continuously  varying  torque  which  is  non-conservative  around  a  cycle  because  of  the 
hysteresis  losses. 

The  digital  computer  program  which  simulates  the  satellite  attitude  dynamics,  and  to  which 
the  magnetic  hysteresis  torque  subroutine  is  added,  is  briefly  described  in  terms  of  its 
flexibilities,  limitations,  and  capabilities  applicable  to  the  present  contract. 

The  various  approaches  used  to  derive  analytical  expressions  to  fit  the  magnetic  curves  are 
described.  Various  methods  had  varying  degrees  of  success.  The  reason  for  abandoning 
each  of  the  unsuccessful  methods  is  given. 

The  work  done  demonstrates  that  the  flux  density  in  a  ferromagnetic  rod  subjected  to  an 
arbitrarily  varying  magnetic  field  can  be  calculated  by  a  digital  computer.  The  limited 
core  storage  and  running  time  economy  dictate  simplifications  in  the  representation  of  the 
magnetic  curves.  A  set  of  ten  polynominals  was  used  to  represent  the  original  data.  In 
the  program,  other  curves  are  interpolated  between  these  polynomials.  An  interpolation 
coefficient  which  varies  with  the  applied  field  is  employed  to  meet  the  various  conditions 
imposed  on  the  interpolated  curves.  From  these  curves,  the  flux  density  in  each  rod  is 
computed  at  every  integration  interval.  These  values  are  then  used  to  compute  the  in¬ 
stantaneous  magnetic  moments  and  hysteresis  torques. 

The  computation  of  other  torques,  calculation  of  vehicle  characteristics  from  input  data, 
integration  of  the  equations  of  attitude  motion  with  respect  to  time,  and  other  required 
computations  are  performed  by  the  existing  digital  program. 
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SECTION  3. 


MAGNETICS  FUNDAMENTALS 

3. 1  FERROMAGNETIC  MATERIALS  AND  HYSTERESIS  LOOPS 

When  a  ferromagnetic  material  is  placed  in  a  magnetic  field,  there  will  be  induced  in  the 
material  a  flux  density  which  is  a  function  of  the  magnitude  of  the  magnetic  field  and  the 
magnetic  characteristics  of  the  material.  This  can  be  expressed  as 


Bi 

=  H.+4irI. 

where 

Bi 

=  internal  magnetic  flux  density 

=  internal  magnetic  intensity 

h 

=  internal  intensity  of  magnetization 

and 

■i 

=  KH. 

where 

K 

=  magnetic  susceptibility. 

(1) 


(2) 


If  the  permeability  of  the  material  is  defined  as 

(3) 

(4) 


u  =1+4  ffK 
=  UH. 


i 


li 


In  a  ferromagnetic  material  the  permeability  is  not  a  constant  but  is  a  function  of  Bi  or 

As  H.  is  varied  from  (H^  max  to  B.  will  vary  from  a  maximum  to  a  minimum  value 

in  a  manner  described  as  a  hysteresis  loop,  such  as  the  outer  loop  in  Figure  3-1.  The 
coordinates  of  the  tips  of  such  a  loop  are  designated  -  h.  (max)  and  t  B(max). 

If  while  traversing  the  hysteresis  loop  from(Hi)max  to  —  (Hj)min,  the  magnitude  of  the 
applied  field  Hj  is  reversed,  a  minor  hysteresis  loop  is  formed.  The  characteristics  of  this 
minor  hysteresis  loop  are  functions  of  the  point  of  flux  reversal.  Further  flux  reversals  will 
produce  additional  minor  hysteresis  loops.  Once  the  major  hysteresis  loop  has  been  traversed 
and  a  field  reversal  takes  place,  the  loop  being  formed  will  bendtowords  the  point  of  last 
reversal.  Internal  reversals  within  a  minor  hysteresis  loop  will  follow  this  property  of 
bending  towards  the  point  of  last  reversal. 

The  magnetic  characteristics  of  any  material  are  a  function  of  the  magnetic  field  intensity 
internal  to  the  material.  In  a  rod,  poles  are  induced  at  the  ends  of  the  rod  which  are 
opposite  to  that  of  the  applied  field.  This  is  a  demagnetizing  field  and  its  effect  is  to  reduce 
the  internal  field  in  the  rod;  it  is  a  function  of  the  length  to  diameter  ratio  of  the  rod. 


H.  -  H0-NI 


(5) 
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where 


Since 


and 


Since 


H.  =  internal  field 
H0  =  external  field 
N  =  demagnetizing  factor 
I  =  intensity  of  magnetization 


Bi  =  H.  +  4  irl, 


I  = 


Bi  -  Hi 
4  v 


H.  =  H  — 
1  o 


H  <<  B. 


~ —  (B.  -  H.) 
4  it  1  » 


(6) 

(7) 

(8) 


N 


H.  =  H  -  , 

l  o  4n  l 


B. 


N 


(9) 


Values  of  are  given  as  a  function  of  and  the  material  permeability  in  Borzorth, 


"Ferromagnetism"  pp.  84G-847.  For  a  long  cylinder  with  infinite  permeability 

4.02  1ogi0(4-)-.92 


(*r): 


(10) 


2  (-7r} 

d 


The  magnetic  moment,  M,  of  a  bar  is 


M  =  ml 

where 

m  =  pole  strength 
l  =  pole  separation 

Since 

3 

II 

It 

AB 

4tt  * 

w  VB 

M  =  — - 
4T7 

(ID 


(12) 


However,  due  to  leakage  and  saturation  effects,  the  pole  separation  is  not  always  equal  to 
the  length  of  the  bar,  Borzorth  and  Chapin*  show  that  the  pole  separation  is  approximately 
.73  the  length  of  the  bar.  Therefore 

M  =  .73  XU  (13) 

4tt 

*Bozorth  and  Chapin,  "Demagnetizing  Factors  of  Rods,"  Journal  of  Applied  Physics, 

Vol.  13,  p.  320,  May,  1942. 
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The  Barkhausen  effect  describes  the  irregularities  in  the  magnetization  of  a  ferromagnetic 
material.  The  magnetization  does  not  take  place  instantaneously  but  in  steps  as  the  field 
strength  is  varied.  As  -gp  approaches  zero,  approaches  the  total  change  in  magnetization 

for  high  nickel  alloys  (~80%  Nj,  20%  Fe).  Investigation  and  consultation  with  Allegheny 

Ludlum  indicates  that  this  should  not  be  a  problem  in  4750  Nickel  Iron. 

It  is  assumed  that  the  rods  are  separated  enc  gh  so  that  there  is  not  magnetic  interaction 
between  them.  Thus  each  rod  acts  as  though  it  is  alone  in  the  magnetic  environment. 

The  rod  material  used  was  AEM  4750  (Allegheny  Electrical  Metal)  which  is  a  nickel  iron 
alloy  containing  approximately  48%  nickel,  the  remainder  being  iron.  In  order  to  obtain 
optimum  magnetic  characteristics  the  material  must  be  annealed.  It  should  be  noted 
that  the  anneal  recommended  by  Allegheny  Ludlum  was  not  that  used  by  the  Applied  Physics 
Laboratory  for  the  rods.  A  comparison  is  given  as  follows: 


Allegheny  Ludlum 

APL 

Heat  Treatment 

2150°F 

1650°F  to  1740°F 

Cooling  Rate 

100°F  per  hour 
max.  to  1100°F 

54°  ±  18°F  per  hour 
from  1100°F  to  570°F 

The  exact  effect  of  this  variation  of  h°at  treatment  is  unknown;  however,  it  is  thought  that 
the  APL  heat  treatment  might  produce  greater  hysteresis  loss  than  the  Allegheny  Ludlum. 

3. 1. 1  Test  Data 

There  is  no  way  to  predict  the  magnetic  properties  of  ferromagnetic  materials  from  any 
fundamental  properties  of  the  material.  The  newest  approach  has  been  the  work  of 
Rayleigh  which  only  applies  to  unsaturated  regions  of  low  flux  density.  The  only  method  for 
obtaining  reliable  data  is  by  test  procedures. 

Allegheny  Ludlum  ran  hysteresis  loops  on  a  ring  sample  of  AEM  4750,  annealed  to  produce 
optimum  magnetic  characteristics;  i.  e. ,  2150°F  anneal.  Curves  were  run  for  Bmax  = 

±  2000,  -  4000,  -  6000  and  -  8000  gauss,  showing  major  and  minor  hysteresis  loops. 

These  curves  are  given  in  Figures  3-1,  3-2,  3-3  and  3-4. 

The  curves  were  obtained  from  a  demagnetized  sample,  first  obtaining  the  major  hysteresis 
loop  and  then  the  minor  hysteresis  loops.  The  slight  shifting  of  the  loops  is  a  result  of 
instability  in  the  test  equipment  rather  than  inherent  magnetic  properties. 

3.  2  CONVERSION  OF  MECHANICAL  ENERGY  TO  HEAT  THROUGH  MAGNETIC 
HYSTERESIS  LOSS 

It  has  been  demonstrated  that  the  torque  acting  on  a  ferromagnetic  rod  placed  in  the  earth's 
magnetic  field  is  associated  with  a  change  in  the  mechanical  energy  equal  and  opposite  to 
the  change  in  the  magnetic  energy  stored  in  the  rod. 


3-3 


-.06 


-.04 


0 


.04 


APPLIED  FIELD  (OERSTEDS) 


Figure  3-1,  Allegheny- Ludlum  2  -  Kilogauss  Data 
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Figure  3-2.  Allegheny- Ludlum  4  -  Kilogauss  Data 
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Figure  3-3.  Allegheny- Ludium  6-Kilogauss  Data 
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Figure  3-4.  Allegheny- Ludium  8-Kiiogauss  Data 
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Figure  3-5  shows  such  a  rod  placed  in  the  earth’s  field,  which  has  the  value  Hg.  The 
component  of  the  field  along  the  rod  is 

Ha  =  Hg  cos  0  (1) 


Figuie  3-5.  Ferromagnetic  Rod  in  the  Earth's  Field 

This  is  the  only  component  which  is  effective  in  producing  any  flux,  <t ,  in  the  rod.  The 
differential  flux  crossing  the  differential  area,  dA,  at  some  point,  is  the  flux  density,  B, 
multiplied  by  the  area, 

d  $  =  BdA,  (2) 

The  flux  density  is 

B  =  u  Hj  =  y(HA  -  DmE),  (3) 

where  Hj  is  the  internal  (induced)  field,  n  is  the  true  permeability  of  the  rod  material, 
and  DM  is  the  demagnetization  factor. 

The  elementary  volume  of  length  d£  is 

dV  =  d  A  d£  (4) 

The  differential  magnetic  moment  of  dV  is 


dM  =  4- 

a 

a 

II 

B  d  A  di 

M 

B  dV. 

(5) 
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The  total  magnetic  moment  of  the  rod  is  the  integral  of  this  expression  throughout  the 
volume,  V, 


M 


=  —  /  BdV. 

4  ir  J 


The  torque  on  the  rod  is 


(6) 


T 


M  Hg  sin  4. 
'if-  HE  sUl  0 


/ 


B  dV. 


(7) 


If  the  rod  rotates  through  an  angle  d9 under  the  influence  of  this  torque,  the  mechanical 
energy  is 


dW  =  T  dfl 
1 


Tv~  HE  sln  e 


d9  | 


B  dV. 


(8) 


The  rotation  will  also  produce  a  differential  change  in  the  applied  field, 


d  =  -Hg  sin  9d  9. 


(9) 


The  corresponding  change  in  the  magnetic  energy,  E,  per  unit  volume,  is 

1 


d2E 


dV 


4  IT 


B  d  Ha. 


(10) 


The  change  in  the  magnetic  energy  of  the  whole  rod  is  the  volume  integral  of  this,  or 


H.  / 


d  E  =  _  d  »A 


Equation  (9)  is  substituted  into  (11) 


B  dV. 


(ID 


d  E  =  -  — 1 —  Hf  sin  9  d  9  J 

4ff  E 


3  dV. 


(12) 


From  equations  (8)  and  (12), 

dW  +  oE  =  0. 


(13) 


This  shows  that  the  mechanical  energy  gained  is  equal  to  the  magnetic  energy  lost,  and 
vice  versa.  The  equation  and  its  method  of  derivation  also  show  that  the  exchange  of 
energy  is  continuous,  because  it  takes  place  with  differential  motions,  It  is  not  necessary 
to  traverse  a  complete  magnetic  hysteresis  loop,  nor  any  specified  fraction  of  it,  for  the 
energy  exchange  to  take  place.  Therefore,  mechanical  energy  is  continually  converted  into 
heat  through  hysteresis  losses  as  the  motion  causes  traversal  of  the  hysteresis  curves. 


3-7 


3.3  COMPUTATION  OF  MAGNETIC  HYSTERESIS  TORQUES 


The  results  of  the  previous  sections  are  used  to  derive  the  torque  expressions  used  in 
the  program  and  the  numerical  value  of  the  conversion  factor. 


Equation  (13)  of  Section  3. 1  relates  the  magnetic  moment  of  one  rod  to  the  flux  density  at 
its  center.  The  factor  .  73  accounts  for  the  decrease  of  flux  toward  the  ends  of  the  rod. 
That  is, 

/ 

B  dV  =  .73  BV.  (1) 


/ 


Each  of  the  rods  used  is  60  inches  long  and  .  11  inches  in  diameter.  The  volume  is  2.  974  n 
cubic  centimeters.  The  ratio  of  the  magnetic  moment  to  the  flux  density  at  the  center  of  the 
rod  is  therefore  .5428  dyne-centimeters  per  oersted-gauss.  This  number  is  doubled  for 
two  rods  (with  negligible  proximity  effect),  and  converted  to  the  units  used  in  the  program. 
TK'  value  is  then  80.067  x  10  foot-pounds  per  oersted-kilogauss.  This  factor  is  designated 
Fj-jC  and  is  an  input  constant  in  the  program.  In  general,  this  factor  accounts  for  the  rod 
length,  rod  diameter,  decrease  of  flux  toward  the  ends  of  the  rod,  number  of  rods, 
proximity  effect,  and  units  conversion  factors. 

The  torque  on  an  elementary  volume,  dV,  of  the  rod  is 

dT  =  I^HE  dV,  (2) 

where  I  is  the  intensity  of  the  magnetization  in  the  rod,  and  H£  is  the  applied  magnetic  field, 
undisturbed  by  the  rod.  The  magnetic  moment  of  the  elementary  volume  is 

dM  =  I  dV.  (3) 


Using  the  expression  for  I  from  equation  (7)  of  Section  3. 1,  with  negligible  in  comparison 
with  P-^,  and  equation  (2)  above,  in  (3),  leads  to 

B. 

dM  =  ^-dV,  (4) 

4  n 

and 

d  T  - 


(5) 


The  torque  on  the  whole  rod  is 

T  =  M  ^he-  (6) 

This  is  consistent  with  equations  (5),  (6),  and  (7)  of  Section  3.2  where  the  value  of  the 
vol  ne  integral  is  given  by  equation  (1)  above. 
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SECTION  4. 

DIGITAL  COMPUTER  SIMULATION 


The  attitude  dynamical  equations  of  the  satellite  were  simulated  by  an  IBM  7094  com¬ 
puter  program  (GOLD-N-ROD),  coded  in  Fortran  II.  This  was  accomplished  by  adding 
a  magnetic  hysteresis  torque  subroutine  to  the  existing  program.  This  program  is 
described  in  detail  in  Reference  1.  The  torque  subroutines  not  needed  for  this  simula¬ 
tion  were  deleted  from  the  program. 

The  applicable  portions  of  the  GOLD-N-ROD  program  are  described  below,  in  sufficient 
detail  to  enable  an  understanding  of  those  capabilities  and  limitations  which  affect  this 
simulation. 

The  magnetic  hysteresis  torque  subroutine  is  described  in  complete  detail.  This  includes 
the  magnetic  curves  for  the  material,  the  fitting  of  analytical  expressions  to  these  curves, 
the  logic  for  tracing  BH  curves  with  an  arbitrarily  varying  H,  and  the  equations  which 
constitute  the  magnetic  hysteresis  torque  subroutine. 

4.  1  GENERAL  DESCRIPTION  OF  THE  DIGITAL  COMPUTER  PROGRAM 

The  program  is  an  IBM  7094  digital  simulation  of  the  attitude  dynamics  of  a  passively 
oriented  satellite.  The  differential  equations  describing  the  attitude  motion  are  numerically 
integrated  with  respect  to  orbit  angle  (in  lieu  of  time).  This  program  was  developed 
prior  -.o  the  NWL  contract.  The  various  features,  capabilities,  and  limitations  of  the 
program  are  best  described  in  connection  v'ith  the  various  subroutines. 

4.1.1  Input  Subroutine  (INPUT) 

For  convenience,  the  various  inputs  are  grouped  in  sets  as  described  below. 

a.  Card  1  -  This  is  the  title  card,  and  may  contain  any  72  alphanumeric  characters. 

b.  Setl.  Earth  Constants  -  At  present,  only  the  solar  pressure  and  solar  heat 
flux  constants  are  listed  as  inputs.  This  is  done  in  order  to  be  able  to  include 
or  exclude  the  effects  of  solar  torque  or  thermal  bendrng  as  desired. 

All  of  the  other  constants  are  internal  to  the  program.  These  include  the 
angle  between  the  ecliptic  and  equatorial  planes,  the  product  of  the  universal 
gravity  constant  and  earth’s  mass,  the  earth’s  orbital  angular  rate,  the  earth's 
spin  rate,  and  the  earth's  equatorial  and  polar  radii.  The  earth’s  orbital 
angular  rate  is  a  constant,  because  the  earth  is  assumed  to  be  in  a  circular 
orbit  about  the  sun.  The  constants  in  the  program  are  in  the  foot-pound-second 
system  and  must  be  changed  if  a  different  system  of  units  is  employed. 

c.  Set  2.  Satellite  Orbital  Parameters  -  These  include  the  geocentric  distances 
at  apogee  and  perigee-  the  precession  rate  and  inclination  of  the  orbit  plane, 
and  the  initial  values  of  the  right  ascension  of  the  ascending  node;  the  time  of 
year,  Greenwich  mean  time;  and  the  orbital  angular  relations  between  the 
ascending  node,  perigee,  and  the  satellite's  position  at  injection. 
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d.  Set  3.  Characteristics  of  the  Central  Body  -  These  include  the  radius  of  the  body 
(assumed  spherical),  its  mass,  moments  of  inertia,  and  the  diffuse  reflectance 
of  its  surface. 

e.  Set?).  Characteristics  of  the  Rod  and  Tip  Mass  -  The  characteristics  of  the 
rod  include  its  length,  radius,  mass  density  per  unit  length,  thickness  oi  wall, 
thermal  conductivity,  linear  coefficient  of  thermal  expansion,  and  the  absorp- 
tance  of  its  surface. 

The  characteristics  of  the  tip  mass  include  its  radius  (the  mass  is  assumed 
to  be  spherical),  the  diffuse  reflectance  of  its  surface,  and  its  mass. 

f.  Set  7.  Initial  Conditions  -  These  are  the  initial  values  of  the  three  Euler  angles 
which  relate  the  body  reference  frame  to  the  orbital  reference  frame,  and  the 
three  body-axis  rates. 

g.  Set 8.  Integration  Parameters,  Time  and  Print  Controls  -  These  include  the 
tolerable  errors  associated  with  the  numerical  integration,  initial  step  size, 
print-out  interval,  and  stop  time.  All  of  these  are  more  fully  explained  in 
the  section  on  numerical  integration. 

4. 1.  2  Initialization  Subroutines  (INITIA  and  VCINIT) 

Certain  parameters  which  are  varied  from  one  run  to  another  are  constant  for  a  parti¬ 
cular  run,  and  are  therefore  initialized  for  economy.  These  include  orbital  parameters, 
elements  of  direction  cosine  matrixes  involving  fixed  angles,  and  vehicle  characteristics. 
The  vehicle  characteristics  include  the  mass,  mass  moments,  and  moments  of  inertia 
of  the  straight  rod  and  certain  combinations  of  parameters  which  are  recurrent  in  the 
solar-torque  equations. 

4.  1.  3  Orbital  Subroutine  (ORBIT) 

The  center  of  mass  of  the  satellite  is  assumed  to  trace  a  circular  or  elliptical  path 
about  the  geocenter.  The  inclination  of  the  orbit  plane  with  respect  to  the  equatorial 
plane  is  constant,  but  provision  is  made  for  a  constant  rate  of  precession  of  the  line  of 
nodes. 

The  outputs  of  this  subroutine  are  time,  satellite  orbital  angle  and  geocentric  distance, 
certain  time  derivatives  of  these  quantities,  altitude  above  the  earth's  surface,  earth's 
orbital  position,  and  a  parameter  used  for  the  earth's  shadow  criterion. 

4.  1.  4  Matrix  Subroutine  (MATRIX) 

This  subroutine  calculates  elements  of  the  direction  cosine  matrixes,  certain  of  their 
time  derivatives,  and  the  satellite’s  latitude  and  longitude. 
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4. 1. 5  Earth’s  Shadow  Criterion 


The  sun's  rays  are  considered  to  be  parallel,  and  therefore  the  earth’s  shadow  is  assumed 
to  be  a  sharply  defined  cylinder.  When  the  satellite  is  in  the  earth’s  shadow,  all  of  the 
solar  effects  are  zero. 

4.1.6  Rod  Geometry  Subroutine  (RODN) 

This  subroutine  is  used  only  when  the  satellite  is  not  in  the  earth’s  shadow.  It  includes 
the  equations  which  describe  the  thermally  bent  rods  under  the  various  conditions  of  full, 
partial,  or  no  illumination.  The  length  <  rod  illuminated  is  dependent  upon  the  shadow 
cast  by  the  (spherical)  central  body,  unless  the  whole  satellite  is  in  the  earth’s  shadow. 
The  mass  moments  of  the  bent  rod  are  computed. 

4.1.7  Center  of  Mass  Subroutine  (CMASS) 

The  instantaneous  moments  and  the  coordinates  of  the  instantaneous  center  of  mass  of 
the  entire  satellite  are  computed.  The  vector  distances  from  this  center  of  mass  to  the 
base  and  tip  of  each  rod  are  computed.  Time  derivatives  of  certain  of  the  foregoing 
variables  are  also  computed. 

4. 1.  8  Moments  of  Inertia  Subroutine  (MOMIN) 

The  instantaneous  moments  and  products  of  inertia  of  the  rod  and  of  the  entire  satellite, 
and  their  time  derivatives,  are  computed  for  the  current  condition  of  rod  illumination. 

4. 1.  9  Solar  Torque  Subroutine  (SOLTOR) 

This  subroutine  is  used  only  when  the  satellite  is  not  in  the  earth's  shadow.  The  solar 
radiation  pressure  torques  on  the  central  body,  rod,  and  tip  mass  are  computed.  The 
sun's  rays  are  all  considered  parallel. 

4. 1. 10  Magnetic  Hysteresis  Torque  Subroutine  (HYSTOR) 

This  subroutine  is  discussed  in  detail  in  section  4.  2. 

4. 1. 11  Gravity  Gradient  Torque  Subroutine  (GRATOR) 

This  subroutine  computes  the  gravity  gradient  torques  on  the  satellite  from  its  moments 
and  products  of  inertia.  The  earth’s  gravitational  field  is  assumed  to  follow  an  inverse 
square  law. 

4. 1. 12  Derivative  Subroutine  (DERIV) 

The  various  torques  are  added  to  obtain  the  total  torques  acting  on  the  satellite.  Euler's 
dynamical  equations  are  solved  explicitly  (in  matrix  form)  to  obtain  the  body-axis  angular 
accelerations.  The  Euler  angular  rates  are  computed  from  the  Euler  angles  and  the 
body-axis  rates.  Each  time  derivate  is  converted  to  the  corresponding  derivative  with 
respect  to  orbital  angle,  by  dividing  by  the  orbital  angular  rate. 
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4. 1. 13  Numerical  Integration  Subroutine  (NOL3) 


The  Adams-Moulton  method  of  numerical  integration  is  used.  A  variable  integration 
interval  is  used.  Specified  error  criteria  are  applied  to  the  estimated  numerical  inte¬ 
gration  error  for  each  variable  to  determine  whether  the  interval  is  satisfactory.  If  the 
interval  is  satisfactory  with  respect  to  every  integrated  variable,  the  results  of  the  inte¬ 
gration  are  accepted  by  the  program.  If  the  interval  is  unsatisfactory,  it  is  halved  until 
a  satisfactory  interval  is  obtained.  Whenever  a  specified  number  of  consecutive  satis¬ 
factory  intervals  ha\  2  been  obtained,  the  interval  is  tentatively  doubled.  The  results 
of  using  the  new  interval  are  tested  against  the  error  criteria  as  explained  above. 

4. 1.  14  Output  S  broutine  (OUTPUT) 

This  subroutine  computes  any  variables  desired  as  output  which  have  not  been  previously 
computed  in  other  subroutines.  Usually,  the  output  time  intervals  are  much  larger  than 
the  integration  intervals,  and  so  the  output  computations  are  performed  only  at  those 
intervals  when  output  is  required.  However,  because  of  the  requirement  for  sometimes 
plotting  the  BH  curves,  more  frequent  outputs  are  desirable  for  the  magnetic  hysteresis 
program.  Therefore,  output  is  obtained  at  every  good  integration  step. 

The  output  variables  required,  which  are  not  available  from  other  subroutines,  include 
the  classical  Euler  angles,  @21  and  63,  one  Euler  angular  rate,  63,  and  the  three 
inertial  coordinates,  X,  Y,  and  Z,  of  the  satellite  position. 

The  Euler  angles  depend  upon  the  values  of  certain  [E]  matrix  elements. 

sin  i>2  =  V^7  •  (1) 

If  this  quantity  is  equal  to  or  ess  than  0.  001,  it  is  considered  to  be  essentially  zero. 

This  yields  the  sing  lar  case,  and  the  computation  and  printout  of  6^,  and  63  are 
omitted  for  that  time  interval.  If  the  quantity  is  greater  than  0.001,  the  computations 
proceed. 

&2  -  cos'1  (En)  .  (2) 

The  computer  subroutine  used  for  the  inverse  cosine  automatically  puts  the  angle  in  the 
first  or  second  quadrant.  is  computed  next. 

if  e13  1  0, 


The  computer  subroutine  used  for  the  inverse  sine  automatically  puts  the  angle  in  the 
first  or  fourth  quadrant. 
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(4) 


II  Ej2  >  0,  and  if  =  0, 

6,  =  180°  -  sin'1  . 

If  Ejg  >  0,  and  if  E^  <  0, 

°1  =  -180°  -  sin’1  (J^-)  . 
0g  is  computed  next. 

if  e31  ^  =  0, 

(^h)  •  . 

If  Egj  <  0,  and  if  E21  -  0, 

6  3  -  180°  *  S1"‘‘  (iST^)  • 
If  Egj  <  0,  and  if  E2j  <  0, 

»3  =  -180°  -  si"'‘ (^V2)  • 
The  angular  rate  0g  is 


(5) 


(6) 


(7) 


(8) 


U.’ 


XI 


i  (E91  a' 


sin 


11  vlJ21  a  Y1  +  ^31  a  Zl 


i )  +  E 


13 


*  > 


(9) 


where  rj  is  the  orbital  angular  velocity,  and  u;^, 
relative  to  an  inertial  frame. 


u'yp  and  are  the  body  axis  rates 


The  inertial  coordinates  of  the  satellite  are  expressed  in  kilometers, 

X  =  -0.3048  x  10"3  Au  Ra,  (10) 

Y  =  -0.3048  x  10'3  A12  Ra,  (11) 

Z  =  0.3048  x  10'3  Al3  Ra,  (1.2) 


where  is  the  geocentric  altitude  of  the  satellite  in  feet,  and  A^,  A^2,  and  A^g  are 
elements  of  the  direction  cosine  matrix  relating  the  orbital  and  inertial  reference  frames. 

4.  2  MAGNETIC  HYSTERESIS  TORQUE  SUBROUTINE  (HYSTOR) 

This  subroutine  computes  the  instantaneous  torque  on  the  satellite  due  to  the  action  of  the 
earth's  magnetic  field  on  the  hysteresis  rods.  e  torque  depends  upon  the  magnetic 
moments  of  the  rods,  which  in  turn  depend  upon  their  magnetic  flux  densities.  The  major 
requirement  of  this  subroutine  is  therefore  the  tracing  of  BH  curves,  or  the  determina¬ 
tion  of  flux  density  under  the  condition  of  a  fluctuating  applied  field.  The  fluctuations 
are  slow,  so  time  rates  of  change  are  negligible. 
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4.  2. 1  Organization  and  Logic  of  'he  Subroutine 


The  subroutine  first  computes  the  components  of  the  earth’s  magnetic  field  in  the  orbital 
(local  vertical)  reference  frame.  This  is  done  by  the  use  jf  a  program  obtained  from  the 
Armed  Forces  Special  Weapons  Command  at  Kirtland  Air  Force  Base,  New  Mexico. 

These  field  components  are  transformed  into  the  satellite  body  reference  frame. 

The  hysteresis  rods  lie  along  the  Y  and  Z  geometric  axes  of  the  satellite  main  body. 

The  same  logic  and  equations  are  used  for  each  of  these  axes.  Therefore,  the  explana¬ 
tion  which  follows  applies  to  each. 

Three  basic  rules  are  used  in  tracing  the  BH  curves,  using  an  arbitrarily  fluctuating 
applied  magnetic  field,  H. 

First,  it  is  assumed  that  the  initial  path  is  the  major  hysteresis  loop  whose  end  points 
correspond  to  the  maximum  values  of  the  earth’s  magnetic  field  encountered  at  the 
perigee  altitude.  Such  a  loop  is  illustrated  in  Figure  4-1,  with  end  points  A  and  B.  The 
assumption  is  justified  on  the  grounds  that  the  satellite  will  undoubtedly  experience  this 
value  of  applied  field  during  its  history  of  assembly,  test,  and  launch. 

Second,  whenever  a  reversal  in  the  algebraic  sign  of  the  field  increment  (equivalent  to 
a  reversal  in  the  sign  of  H)  occurs,  the  new  BH  curve  must  pass  through  the  last  prior 
reversal  point.  In  Figure  4-1,  a  reversal  occurs  at  point  C,  and  the  curve  CD,  if  extended, 
would  pass  through  point  B.  Similarly,  the  curve  DE,  if  extended,  would  pass  tnrough 
point  C,  etc.  Thus  the  end  point  of  any  curve  is  the  beginning  point  of  the  previous  curve. 
This  rule  is  based  on  the  known  effect  of  applying  an  alternating  field  superimposed  on  a 
direct  field.  That  is,  a  minor  loop  between  the  same  two  end  points  is  traced  repetitively. 
In  the  case  of  a  rapidly  alternating  field,  several  cycles  might  be  required  before  the 
loop  becomes  evactly  repetitive.  However,  it  is  believed  that  with  slow  variations  in  the 
applied  field,  the  first  traverse  of  the  loop  very  closely  approximates  what  would  become 
the  final  path.  It  should  be  noted  that  there  is  no  implication  here  that  a  complete  minor 
loop  is  actually  traced,  but  only  that  the  path  of  the  BH  point  is  along  a  curve  which 
passes  through  the  prior  reversal  point.  Of  course,  if  the  variation  in  H  is  great  enough, 
that  point  will  be  reached,  in  which  case  the  complete  minor  loop  will  have  been  traced 
once. 

Third,  whenever  the  value  of  H  ranges  beyond  the  end  point  of  the  present  curve  (the 
prior  reversal  point  discussed  above),  the  path  reverts  to  the  previous  curve  of  like  kind 
(that  is  the  previous  ascending  or  descending  curve).  For  example,  if  the  BH  point  is  on 
curve  FE  of  Figure  4-1,  and  H  subsequently  increases  beyond  E,  the  curve  DC  will  then  be 
followed.  This  rule  may  not  be  strictly  correct,  but  it  is  the  concensus  of  expert  opinion 
that  it  is  a  good  approximation,  and  as  good  as  can  be  done  with  available  knowledge. 
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The  implementation  of  these  basic  rules  is  illustrated  by  the  logic  chart  of  Figure  4-2. 

The  two  disconnected  portions  of  this  chart  are  a  result  of  the  nature  oi  the  numerical 
integration  subroutine  used  in  the  GOLD-N-ROD  program.  This  subroutine  achieves  a 
balance  between  accuracy  and  economy  by  choosing  the  integration  interval  so  that  the 
estimated  errors  meet  a  set  of  specified  erroi  criteria.  Only  if  the  criteria  are  satisfied 
is  the  interval  considered  a  good  integration  step.  It  is  only  at  such  ?.  step  that  the  logic 
based  on  the  above  rules  is  applied.  However,  the  integration  method  requires  tne  com¬ 
putation  of  derivatives,  and  therefore  of  torques,  at  every  interval  and  at  certain  sub¬ 
intervals.  At  ever,,  interval  or  subinterval,  tne  calculation  of  magnetic  moments  and 
torques  is  based  on  the  value  of  the  flux  density  at  the  most  recent  good  step. 

At  the  initial  time,  the  BH  point  is  assumed  to  be  on  the  major  hysteresis  loop,  ar 
stated  previously.  In  order  to  determine  whether  the  point  lies  on  the  ascending  01 
descending  branch  of  this  loop,  the  initial  values  of  the  instantaneous  time  derivatives 
of  the  earth's  field  components  are  computed.  If  the  derivative  in  question  has  a  positive 
value,  the  initial  point  lies  on  the  ascending  curve  and  the  ascending-descending  index 
(D-A  index)  is  set  equal  to  +1.  Conversely,  if  the  initial  value  of  H  is  negative,  the 
initial  point  is  on  the  descending  branch,  and  the  D-A  index  is  set  equal  to  -1. 

In  order  to  use  the  reversal  criterion  (as  explained  later)  at  the  next  good  integration 
step,  an  artificial  increment  in  the  field  must  be  created,  which  corresponds  to  the  time 
increment  prior  to  the  initial  time.  Only  the  algebraic  sign  of  this  artificial  increment 
is  used  in  the  logic,  and  the  numerical  value  is  of  no  consequence. 

Each  BH  curve  is  identified  by  the  coordinates  of  its  beginning  point.  The  end  point  of 
the  curve  is  the  beginning  point  of  the  prior  curve,  in  accordance  with  the  rules  previously 
explained.  Because  of  the  necessity  of  sometimes  returning  to  a  prior  curve,  the  be¬ 
ginning  point  coordinates  of  every  curve  must  be  stored,  unless  they  are  intentionally 
erased  at  a  later  time.  The  section  of  the  memory  serving  this  function  is  called  the 
curve  array.  The  beginning  point  coordinates  of  each  curve  are  identified  by  the  curve 
array  index  (JHYS),  which  is  assigned  in  chronological  order.  In  accordance  with  this 
system,  the  branch  of  the  major  loop  on  which  the  initial  point  lies  is  assigned  JHYS 
equal  to  2,  and  the  beginning  point  coordinates  are  stored  in  the  curve  array  with  this 
index  value.  The  other  branch  of  the  major  loop  is  assigned  JHYS  equal  to  1,  and  its 
beginning  point  coordinates  are  stored  accordingly.  Henceforth,  the  beginning  point  of 
the  curve  presently  used  always  has  an  index  equal  to  JHYS,  the  highest  value  of  the 
index.  The  end  point  of  the  present  curve,  being  the  beginning  point  of  the  prior  curve, 
has  an  index  equal  to  JHYS-1. 

At  the  initial  time,  the  flux  density  in  the  rods  is  computed  from  the  magnetic  field  com¬ 
ponents  (in  the  satellite  body  frame),  using  the  analytical  approximation  for  the  major 
loop,  which  has  been  stored  as  one  of  the  eleven  reference  curves. 

At  any  interval,  including  the  initial  time,  the  magnetic  moments  of  the  rods  are  com¬ 
puted  from  their  flux  densities.  After  th's  has  been  done  for  both  the  Y  and  Z-axis  rods, 
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Figure  4-2.  BH  Curve  Tracing  Logic 
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the  values  are  added  to  the  respective  components  of  the  fixed  or  residual  magnetic 
moment  of  the  satellite  to  obtain  the  components  of  the  total  magnetic  moment.  (The 
fixed  magnetic  moments  of  the  NWL  satellite  are  negligible. )  The  components  of  the 
magnetic  torque  are  then  computed  from  those  of  the  satellite  magnetic  moment  and 
those  of  the  earth's  magnetic  field. 

At  any  good  integration  step,  it  is  necessary  to  detect  a  reversal  in  the  sign  of  the  in¬ 
crement  of  the  field,  H.  Such  reversals  occur  at  points  C,  D,  E,  and  F,  in  Figure  4-1. 
Detection  of  reversals  is  accomplished  by  comparing  the  last  two  increments.  If  they 
are  of  like  sign,  as  manifested  by  a  positive  product,  no  reversal  has  occurred.  With 
no  reversal,  the  D-A  index  remains  the  same  as  previously.  If  the  last  two  increments 
in  H  are  of  unlike  sign,  as  manifested  by  a  negative  product,  a  reversal  has  occurred. 

The  sign  of  the  D-A  index  is  then  reversed  and  checked.  A  reversal  also  has  occurred 
if  the  last  increment  in  H  is  zero,  as  manifested  by  a  zero  product.  The  procedure 
illustrated  in  the  chart  treats  this  case  as  though  the  product  was  negative.  This  pre¬ 
vents  the  occurrence  of  a  zero  product  at  the  next  interval,  which  would  otherwise  cause 
a  second  apparent  reversal  for  a  single  occurrence  of  a  zero  increment.  The  logic  shown 
in  the  chart  will  correctly  determine  the  occurrence  of  reversals  at  the  following  time 
interval  for  any  sign  (or  zero  value)  of  the  increment  pertaining  to  that  interval. 

After  any  reversal,  a  new  BH  curve  must  be  derived,  and  is  then  treated  as  the  present 
curve.  This  is  done  by  using  tne  coordinates  of  the  beginning  and  end  points,  as  explained 
later.  The  beginning  point  coordinates  are  the  values  of  B  and  H  at  the  time  interval 
prior  to  the  one  when  the  reversal  was  detected.  These  prior  values  have  always  been 
stored,  to  facilitate  this,  because  it  is  never  known  when  the  next  interval  will  produce 
a  reversal.  These  beginning  point  coordinates  are  assigned  the  next  available  JHYS 
number  and  stored  in  the  curve  array.  The  coordinates  of  the  end  point  of  the  new  <  urve 
are  the  values  then  available  in  the  curve  array  under  the  index  JHYS-1. 

Whether  a  reversal  has  occurred  or  not,  the  value  of  H  must  be  tested  to  determine 
whether  it  is  beyond  the  end  point  of  the  present  BH  curve.  If  not,  and  if  no  reversal  has 
occurred,  the  flux  density  is  computed  by  using  the  present  curve.  If  a  reversal  has  just 
occurred,  but  H  is  not  beyond  the  end  point  of  the  new  curve,  that  new  curve  is  used  in 
determining  the  flux  density. 

Whenever  the  value  of  H  lies  beyond  the  end  point  of  the  present  curve,  regardless  of 
the  occurrence  of  a  reversal,  the  last  two  curves  must  be  erased  from  the  memory 
(curve  array),  and  the  curve  prior  to  those  two  must  be  considered  for  use.  Since  it  is 
possible  that  the  value  of  H  lies  beyond  the  end  point  of  any  new  curve,  the  test  must 
be  repeated,  as  shown  in  the  chart,  until  a  satisfactory  end  point  has  been  determined. 

The  end  point  is  accepts  >le  if  the  value  jf  H  does  not  lie  beyond  it.  The  coordinates  of 
the  acceptable  end  point,  together  with  those  of  the  beginning  point  of  the  corresponding 
curve,  are  used  to  derive  a  new  BH  curve  (which  then  becomes  the  present  curve)  by  a 
method  explained  later.  The  flux  density  is  then  computed  from  this  new  curve. 
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In  the  program,  all  B-H  curves  are  considered  as  ascending  curves.  When  the  path  of 
the  BH  point  is  along  a  descerding  curve,  it  is  considered  for  computational  purposes 
that  the  signs  of  B  and  H  are  reversed  for  the  entire  history  of  the  rods.  Then  the 
descending  curve  becomes  an  ascending  one.  The  computed  value  of  the  flux  density 
is  multiplied  by  the  D-A  index.  The  use  of  this  index  in  appropriate  places  automatically 
inverts  curves  and  multiplies  computed  values  by  -1  whenever  necessary,  in  order  to 
use  the  basic  ascending  curves  and  associated  logic  for  both  ascending  and  descending 
curves. 

At  every  good  integration  step,  the  variables  required  for  printout  are  stored  on  tape. 

4.  2.  2  Inputs 

The  inputs  for  the  hysteresis  torque  subroutine  include  two  constants  and  88  coefficients. 
One  constant  is  the  value  of  the  external  magnetic  field  corresponding  to  the  tip  of  the 
major  hysteresis  loop.  The  second  is  a  conversion  factor  to  compute  the  magnetic 
moments  from  the  flux  densities.  The  88  coefficients  are  those  of  eleven  sev  ith-degree 
polynomials,  which  are  used  to  represent  the  BH  curves.  These  are  discussed  in  detail 
in  Section  4.  3. 1. 

4.  2.  3  Equations 

Certain  of  the  computations  are  performed  at  every  integration  step  and  others  only  at 
good  integration  steps.  This  is  accomplished  by  entering  the  HYSTOR  subroutine  from 
other  subroutines.  At  every  integration  step,  HYSTOR  is  entered  from  DERIV,  with  the 
entry  point  index,  NOPT,  equal  to  1.  The  computations  described  in  Section  4.  2.  3. 1 
are  then  performed.  At  good  steps  only,  HYSTOR  is  entered  from  OUTPUT,  with  the 
entry  point  index  equal  to  2.  The  computations  described  in  Section  4.  2.  3.  2  are  then 
performed. 

4.  2.  3. 1  Computations  at  Every  Integration  Step 

At  every  integration  step,  certain  of  the  following  computations  are  performed.  At 
every  time  interval,  the  components  of  the  earth's  magnetic  field  are  computed,  as 
described  in  Section  4,  2.  3.  la.  At  the  initial  time,  TQ,  the  computations  described  in 
Section  4.  2.  3.  lb  are  done.  At  every  time  interval,  the  magnetic  moments  and  torques 
are  computed,  as  described  in  Section  4.  2.  3.  lc. 

The  time,  T,  is  tested.  If 

T  -  Tq  -  0. 01  5  0,  •  (1) 

the  time,  T,  is  considered  equal  to  the  initial  time,  T  .  The  small  increment  assures 
that  the  correct  choice  will  be  made,  regardless  of  roundoff  errors. 
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a.  Components  of  the  Earth's  Magnetic  Field  -  The  components  of  the  earth’s 
magnetic  field,  in  spherical  coordinates,  B^,  B^,  and  BM,  are  computed  in 
a  subroutine  adapted  from  the  program  used  by  the  Armed  Forces  Special 
Weapons  Command  at  Kirtland  Air  Force  Base,  New  Mexico.  This  subroutine 
requires  the  altitude,  A^,  latitude,  a,  and  longitude,  Qy,  of  the  satellite, 
as  inputs.  These  are  all  available  from  previous  subroutines.  The  components 
in  spherical  coordinates  are  converted  to  those  in  the  orbital  and  satellite 
reference  frames.  Use  is  made  of  the  right  ascension  of  the  satellite,  QT, 
and  the  right  ascension  of  the  ascending  node,  also  available  from  previous 
subroutines.  These  celestial  angles  are  measured  from  a  baseline  along  the 
vernal  equinox.  The  positive  direction  of  this  baseline,  corresponding  to  zero 
values  of  the  angles,  is  from  the  heliocenter  to  the  geocenter  at  the  time  of 
vernal  equinox.  This  is  opposite  from  the  customary  reference  direction. 

The  difference  in  the  longitudes  of  the  satellite  and  of  the  ascending  node  is 

^M  =  ^T  ”  ^N’  (2) 


The  [D]  matrix  is  used  to  convert  from  the  spherical  coordinate  frame  to 
the  orbital  frame,  RPQ.  The  elements  of  the  [D]  matrix  are 


D11 = 

(3) 

D22  =  sin  v  cos  nM, 

(4) 

t'v  COS  V 

23  ‘  ‘  cos  x  ’ 

(5) 

D32  =  "D23» 

(6) 

D33  =  °22’ 

(7) 

D12  =  Dl3  =  D21  =  °31  =  °» 

(8) 

where  v  is  the  orbital  inclination  angle. 

The  components  of  the  earth’s  magnetic  field,  in  the  orbital  frame,  are 


'hr‘ 

■-Br 

HP 

II 

f - 1 

o 

1 _ J 

B6 

.hqJ 

B. 

The  [E]  matrix  is  used  to  convert  from  the  orbital  frame  to  the  satellite  frame. 
The  elements  of  this  matrix  are  computed  in  a  previous  subroutine.  The  com¬ 
ponents  of  the  field,  in  the  satellite  frame,  are 
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(10) 


’H!3_ 

V 

H11 

-  [E] 

HP 

K12 

hq 

is  the  component,  is  the  Y^  component,  and  is  the  Zj  com¬ 
ponent. 


b.  Computations  at  the  Initial  Time,  Tp  -  The  logic  requires  the  time  increment 
in  the  magnetic  field.  At  the  initial  time,  T  ,  the  field  increment  would 
correspond  to  the  time  interval  prior  to  T  .  This  increment  is  not  available, 
and  therefore  it  must  be  generated  artificially.  Only  the  algebraic  sign  of  the 
increment  has  any  effect,  not  the  numerical  value. 

The  increments  due  to  satellite  attitude  motion,  satellite  orbital  motion,  and 
earth  spin  motion  are  computed  and  combined.  The  increments  due  to  attitude 
motion  depend  on  the  time  derivatives  of  the  [E]  matrix  elements.  These  are 
designated  by  the  additional  subscript,  D: 


D21  = 

E31 

^Xl  ' 

■Ell 

ojZi  + 

E22 

b  , 

(ID 

D22  = 

E32 

“xi  ‘ 

'  E12 

a;Zl  ' 

E21 

V  , 

(12) 

D23  = 

E33 

aXl  ' 

‘  El3 

aZl’ 

(13) 

D31  = 

E11 

^  Y1  ' 

‘  E21 

a  XI  + 

E32 

V  , 

(14) 

D32  = 

E12 

o-'yi  - 

'  E22 

aXl  " 

E31 

*  . 

(15) 

D33  = 

E13 

wyi  • 

E23 

“xp 

(16) 

where  r?  is  the  orbital  angular  velocity,  and  the  body  axis  rates  are  u.^, 
u*Yp  and  u.1  ^ 

The  required  time  derivatives  of  the  components  of  the  earth's  magnetic  field 
are 


H  D1  "  ED21  HR  +  ED22  HP  +  ED23  HQ’ 

H1D2  =  ED3l  HR  +  ED32HP  +  ED33  HQ’ 
or 

H1DI  =  ED,I+1,1HR+  ED,  1+1,2  HP  +  ED,  1+1,3  HQ- 

The  subscript  I  is  1  for  the  Y  axis  and  2  for  the  Z  axis.  The  same  set  of 
equations  is  used  for  both  axes,  in  most  cases.  Only  one  set  of  equations, 
with  subscript  I,  will  be  included  hereafter. 


(17) 

(18) 
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The  increments  in  the  field  due  to  the  satellite  orbital  motion  and  earth  spin 
motion  depend  upon  extrapolated  values  of  the  satellite's  latitude  and  longitude 
prior  to  T  .  The  Tq  value  of  the  altitude  is  used,  because  the  change  would 
be  so  small,  anyway. 

The  orbital  angle  prior  to  time  zero  is 

T7d  =  r?  -  ,  (20) 

where  Ap  is  the  initial  integration  interval.  The  extrapolated  prior  latitude 
is 

A  j-j  =  sin"1  (sin  ^  sin  r?D).  (21) 


The  corresponding  longitude  is  computed  from  various  angles.  The  sidereal 
angle  of  Greenwich  is 

nGD  =  QG  "  ^G  ^ 

TJ 

where  CiG  is  the  value  at  T  f>Q  is  the  earth's  spin  rate,  and  n  is  the  satellite's 
orbital  angular  velocity. 


The  right  ascension  of  the  ascending  node  is 


QND  =  nN  ’  nN 


An 

n 


(23) 


where  oN  is  the  value  at  Tq  and  is  the  precession  rate  of  the  nodes.  The 
rignt  ascension  of  the  satellite,  is  found  from 


sin  O  ipj-j  - 


sin  Ojsjj)  cos  +  cos  v  cos  sin 
cosT^ 


(24) 


cos  Cl 'pjj  ■ 


cos  Qjsjq  cos  -  cos  v  sin  0^  sin 
cos  A.p 


(25) 


The  extrapolated  prior  longitude  is 
Gvd  =  ftTD  ‘  fiGD  • 


(26) 


The  corresponding  difference  in  the  longitudes  of  the  satellite  c.nd  of  the 
ascending  node  is 


nMD  “  RTD  *  °ND  • 


(27) 
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These  values  are  used  in  the  earth's  magnetic  field  subroutine  to  find  the 
extrapolated  prior  values  of  the  field  in  spherical  coordinates.  The  elements 
of  the  [D]  matrix  are  found  from  Equations  (3)  through  (8),  with  used 

in  place  of  0^,  and  X^  used  in  place  of  x.  Equations  (9)  and  (10)  are  then 
ufed  to  compute  Dp^g  ^EHll’  anc*  ^EH12’  the  comPonents  °f  the  extrapolated 
prior  field  in  the  satellite  reference  frame.  The  artificial  increments,  due  to 
orbital  travel  and  earth’s  spin,  are 

DH11  =  H11  ‘  DEH11’ 

DH12  =  H12  "  DEH12  ‘  (29* 

The  total  artificial  increments  are 

HHNI  =  H1DI  '  DH1I  (30) 

The  D-A  index,  is  determined  in  accordance  with  the  algebraic  sign  of 

HHNI 

U  HHNI  °»  HDAI  =  'L 

H  HHNI  -  °>  HDAI  =  (32) 

The  end  point  coordinates  of  the  branch  of  the  major  hysteresis  loop  are 

HHC1I  "  HDAI  HFAH’  (33) 

HBC1I  =  "HDAI  HFA11»  (34) 

where  Hp^  and -Hp^j  are  the  coordinates  of  the  positive  tip  of  the  loop. 

The  corresponding  beginning  point  coordinates  are 

HHC2I  =  "HHC1I»  (35) 

HBC2I  =  'hbcii-  (36) 

The  argument  of  the  polynomial  used  to  find  the  flux  density  is  the  D-A -indexed 
difference  between  the  field  value  and  the  beginning  point  value, 

X  =  HDAI  (H1I  "  HHC2I)*  (37) 

The  value,  Hp,  of  the  polynomial  is  computed.  The  flux  density  is 


HBLI  =  hdai  hf- 
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The  increments  in  the  field  due  to  the  satellite  orbital  motion  and  earth  spin 
motion  depend  upon  extrapolated  values  of  the  satellite’s  latitude  and  longitude 
prior  to  T  .  The  TQ  value  of  the  altitude  is  used,  because  the  change  would 
be  so  small,  anyway. 

The  orbital  angle  prior  to  time  zero  is 

=  V  -  ar?  ,  (20) 


where  At?  is  the  initial  integration  interval, 
is 

Ap  -  sin”1  (sin  sin  r?D). 


The  extrapolated  prior  latitude 

(21) 


The  corresponding  longitude  is  computed  from  various  angles.  The  sidereal 
angle  of  Greenwich  is 


*GD  “  aG  ” 


AIL 

n 


(22) 


where  fiG  is  the  value  at  Tq,  qg  is  the  earth’s  spin  rate,  and  tj  is  the  satellite's 
orbital  angular  velocity. 


The  right  ascension  of  the  ascending  node  is 


°ND  =  nN  “  aN 


At? 

v 


(23) 


where  nN  i*>  the  value  at  Tq  and  nN  is  the  precession  rate  of  the  nodes.  The 
right  ascension  of  the  satellite.  is  found  from 


sin  Q  -pj-j  - 


sin  nND  cos  t?d  +  cos  v  cos  sin 
_  -  cos  X  D 


(24) 


cos  0TD  = 


C°S  qnd  cos  r?D  -  cos  V  sin  0ND  sin  r?D 
cos  A  D 


(25) 


The  extrapolated  prior  longitude  is 
nVI)  =  ftTD  ”  °GD  * 


(26) 


The  corresponding  difference  in  the  longitudes  of  the  satellite  and  of  the 
ascending  node  is 

nMD  "  ^TD  ”  °ND  ' 
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These  values  are  used  in  the  earth's  magnetic  field  subroutine  to  find  the 
extrapolated  prior  values  of  the  field  in  spherical  coordinates.  The  elements 
of  the  [D]  matrix  are  found  from  Equations  (3)  through  (8),  with  used 

in  place  of  0^,,  and  Xq  used  in  place  of  a.  Equations  (9)  and  (10)  are  then 
used  to  compute  ^EHll’  ancl  ^EH12  the  comPonents  of  the  extrapolated 

prior  field  in  the  satellite  reference  frame.  The  artificial  increments,  due  to 
orbital  travel  arid  earth's  spin,  are 

°H11  =  H11  ‘  DEH11’  (28) 

DH12  =  H12  '  DEH12  •  <29) 

The  total  artificial  increments  are 

HHNI  =  H1DI  *  °H1I  (30) 

The  D-A  index,  H^j  is  determined  in  accordance  with  the  algebraic  sign  of 
HHNI 

U  HHNI  °’HDAI  =  _1, 

U  HHN1  ?  °»  HDAI  =  +1*  (32) 

The  end  point  coordinates  cf  the  branch  of  the  major  hysteresis  loop  are 

HHC1I=  HDAI  HFAH’  (33) 

HBC1I  =  "HDAI  HFA1P  (34) 

where  anc*  "^fAII  are  the  coordinates  of  the  positive  tip  of  the  loop. 

The  corresponding  beginning  point  coordinates  are 

HHC21  =  “HHC1I»  (35) 

HBC2I  =  _HBCir  (36) 

The  argument  of  the  polynomial  used  to  find  the  flux  density  is  the  D-A -indexed 
difference  between  the  field  value  and  the  beginning  point  value, 

X  =  HDAI  (H1I  ‘  HHC2I)'  (37) 

The  value,  Hp,  of  the  polynomial  is  computed.  The  flux  density  is 


HBLI  "  hdai  hf- 
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The  increments  in  the  field  due  to  the  satellite  orbital  motion  and  earth  spin 
motion  depend  upon  extrapolated  values  of  the  satellite’s  latitude  aim  longitude 
prior  to  T  .  The  Tq  value  of  the  altitude  is  used,  because  the  change  would 
be  so  small,  anyway. 

The  orbital  angle  prior  '.o  time  zero  is 

bD  =  n  -  *v  ,  (20) 


where  Arj  is  the  initial  integration  interval.  The  extrapolated  prior  latitude 
is 

a  p  =  sin-1  (sin  1/  sin  rjp).  (21) 


The  corresponding  longitude  is  computed  from  various  angles.  The  sidereal 
angle  of  Greenwich  is 


nGD  "  °"G  ‘  ^G 


AH 

v 


(22) 


where  is  the  vaiue  at  Tq,  Qq  is  the  earth’s  spin  rate,  and  i?  is  the  satellite’s 
orbital  angular  velocity. 

The  right  ascension  of  the  ascending  node  is 

nND  =  aN  '  nN  ’  (23) 


where  nN  is  the  value  at  Tq  and  is  the  precession  rate  of  the  nodes.  The 
right  ascension  of  the  satellite,  is  found  from 


sin  n  "p£)  - 


sin  cos  +  cos  u  cos  (1^^  sin 

COS  *D 


(24) 


COS  fl 'pjj  - 


cos  cos  Tj -  cos  v  sin  sin  r?p 
cos  A  p 


(25) 


The  extrapolated  prior  longitude  is 
nVD  =  fiTD  ‘  fiGD  ’ 


(26) 


The  corresponding  difference  in  the  longitudes  of  the  satellite  and  of  the 
ascending  node  is 


nMD  “  ^TD  ‘  °ND  ’ 


(27) 
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These  values  are  used  in  the  earth’s  magnetic  field  subroutine  to  find  the 
extrapolated  prior  values  of  the  field  in  spherical  coordinates.  The  elements 
of  the  [D]  matrix  are  found  from  Equations  (3)  through  (8).  with  used 

in  place  of  0^,  and  A  ^  used  in  place  of  a.  Equations  (9)  and  (10)  are  then 
used  to  compute  Dp,^  ^EHll’  anc*  ^EH12’  the  comPonents  the  extrapolated 
prior  field  in  the  satellite  reference  frame.  The  artificial  increments,  due  to 
orbital  travel  and  earth's  spin,  are 

DH11  =  H11  "  DEH11’  *28) 

DH12  =  H12  "  DEH12  *  ^29) 

The  total  artificial  increments  are 

HHNI  =  H1DI  ^  +  °H1I  (30) 

The  D-A  index,  H^j  is  determined  in  accordance  with  the  algebraic  sign  of 
HHNI 

K  HHNI  °»  HDAI  =  ~1‘  *3l) 

U  HHNI  "  °»  HDAI  =  +1*  (32) 

The  end  point  coordinates  of  the  branch  of  the  major  hysteresis  loop  are 

HHC1I  =  HDAI  HFAH’  (33) 

hbcii  =  -hdai  HFA11»  *34) 

where  Hp^H  and are  the  coordinates  of  the  positive  tip  of  the  loop. 

The  corresponding  beginning  point  coordinates  are 

HHC2I  =  _HHC1I»  (35) 

HBC2I  =  "hbci:*  ^36) 

The  argument  of  the  polynomial  used  to  find  the  flux  density  is  the  D-A-indexed 
difference  between  the  field  value  and  the  beginning  point  value, 

X  =  HDAI  (H1I  “  HHC2I)*  ^37) 

The  value,  Hp,  of  the  polynomial  is  computed.  The  flux  density  is 


HBLI  "  HDAI  hf* 
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The  vaiue,  Hjj,  of  the  magnetic  field  component  is  assigned  to  the 

symbol  for  the  prior  value  so  that  it  will  be  available  as  such  at  the  next 
time  interval. 


HHMI  ~  Hir 


(39) 


c.  Magnetic  Moments  and  Torques  -  The  magnetic  moments  and  torques  are 
computed  at  every  time  interval.  The  magnetic  moments  are  computed  from 

HDBY  =  HDCHBL1  *  HCY’  (40) 

HDBZ  =  HDCHBL2  +  HCZ’  <41) 

where  Hj^  is  the  constant  for  converting  from  flux  densiiy  at  the  center  of 
one  rod  to  total  magnetic  moment  of  a  set  of  two  parallel  rods.  For  the  NWL 

_r* 

satellite,  this  factor  is  80.  067  x  10  foot-pounds  per  oersted-gauss.  H^, 
H^y,  and  H^z  are  the  components  of  the  fixed  or  residual  magnetic  moment 
of  the  satellite.  For  the  NWL  satellite,  this  magnetic  moment  is  negligible. 

The  components  of  the  magnetic  hysteresis  torque  are 


TMHX  =  HDBYH12  ‘  HDBZH11’ 
TMHY  =  HDBZ  H13  "  HCX  H12’ 
TMHZ  =  HCXH11  '  HDBYH13- 


(42) 

(43) 

(44) 


These  torques  are  added  to  the  other  torques  acting  on  the  satellite. 
The  program  then  returns  to  the  DERIV  subroutine. 


4.  2.  3.  2  Computations  at  Good  Integration  Steps 

At  good  integration  steps  only,  the  logic  (which  implements  the  basic  rules  for  tracing 
hysteresis  loops)  is  employed. 

The  time,  T,  is  tested  by  inequality  (1).  If  this  inequality  holds,  T  is  considered  equal 
to  Tq,  and  the  curve  array  indexes  for  both  axes  are  set  equal  to  2, 


JHYS1  "  2’ 

(45) 

J  _  9 

HYS2  ‘ 

(46) 

The  program  then  returns  to  the  output  subroutine,  and  the  required  variables  are 
stored  on  tape  for  later  printout. 
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If 


T  -  Tq  -  0.  01  >  0, 


(47) 


the  time  is  later  than  the  initial  time,  and  the  ensuing  logic  is  employed  for  both  axes 
by  the  index  I  taking  the  values  1  and  2. 

The  index  is  set  equal  to  zero  for  possible  later  use. 

The  new  value  of  the  previous  field  increment  is  set  equal  to  the  old  value  of  the  present 
increment 


HHPI  “  HHNT 


(48) 


The  new  value  of  the  present  increment  is  computed  from  the  appropriate  values  of  the 
field. 


HHNI  =  H1I  '  hhmi- 


(49) 


Writing  the  new  values  over  the  old  ones  automatically  erases  the  latter  from  the  memory. 


In  a  similar  manner,  the  new  value  >  i  the  previous  flux  density  is  set  equal  to  the  old 
present  value, 


HBMI  =  hbli- 


(50) 


The  product  of  the  present  and  previous  field  increments  is  computed, 


HHQI  '  HHNIHHPF  (51) 

If  this  product  is  positive,  no  reversal  in  H  has  occurred,  and  the  procedure  of  Section  b 
below  is  followed. 

If  the  product  is  negative,  a  reversal  in  H  has  occurred,  and  the  procedure  of  Section  a 
below  is  followed. 

If  the  product  is  zero,  a  reversal  in  H  has  occurred.  In  this  case  the  present  increment 
is  set  equal  to  the  negative  of  the  previous  increment, 


HHNI  ~  'hhpi- 


(52) 


and  the  procedure  of  Section  a  below  is  followed. 

a.  Procedure  With  a  Reversal  of  H  -  When  a  reversal  of  H  occurs,  the  sign  of 
the  D-A  index  is  reversed, 


HDAI  "  *HDAI’ 


(53) 


The  new  value  of  the  index  is  checked  for  agreement  with  the  sign  of  *he 
present  field  increment.  If 


HDAI  HHNI  °» 


(54) 


i 
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the  program  is  stopped.  If 


HDAI  HHNI  ~  °’ 


(55) 


the  procedure  continues. 

The  curve  array  index  is  checked  next.  If 

JHYSI  ‘  400  =  °<  (56) 

the  remaining  storage  is  insufficient,  and  the  program  is  stopped.  If 

JHYSI  '  400  ^  °>  (57) 


sufficient  storage  remains,  and  the  program  continues. 
The  curve  array  index  is  increased  by  one, 

J  =  JHYSI  +  l» 

JHYSI  =  J* 


(58) 

(59) 


The  new  beginning  point  coordinates  are  set  equal  to  the  previous  values  of 
the  field  and  flux  density. 


HHCJI  "  HHMI' 
HBCJI  =  HBMI 


(60) 

(61) 


These  coordinates  are  automatically  stored  in  the  curve  array. 


The  index  is  set  equal  to  1,  as  an  indication  that  a  new  BH  curve  must 

be  used. 


The  value  of  the  field  is  tested  to  determine  whether  it  exceeds  the  new  end 
point.  To  facilitate  this,  a  temporary  index  is  needed, 

•J  =  jhySF  (62) 


The  H -coordinate  of  the  new  end  point  is  then  j  i  i 
difference  between  the  field  and  this  coordinate  is 


The  D-A -indexed 


TEMPI  -  HDAI  (H1I  "  HHC,  J-1,1** 


(63) 


If 

TEMPI  =  °’ 


(64) 
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the  field  exceeds  the  new  end  point.  The  curve  array  index  is  then  reduced 
by  two, 


JHYSI  ~  JHYSI  "  2* 


(65) 


The  index  N^gg-pj  is  again  set  equal  to  1.  Since  J^ySI  was  rec*uced,  it  is 
advisable  to  check  that  its  value  is  not  too  low.  The  minimum  allowable  value 
is  two.  If 

•IrysI  '  ^  (66) 

the  index  is  one  or  less,  and  the  program  is  stopped.  If 

Juvci  “1^0,  (67) 


the  index  is  at  least  two,  and  the  program  continues. 

The  reduction  in  the  index  results  in  the  admission  of  a  new'  BH  curve,  with 
new'  beginning  and  end  points.  The  test  of  the  field  for  exceeding  the  new  end 
point  is  repeated  by  again  using  Equations  (62)  and  (63).  As  long  as  inequality 
(64)  holds,  the  repetition  is  continued.  Eventually,  when 

TEMPI  '  °»  (68) 


a  curve  has  been  identified  whose  end  point  is  not  exceeded  by  the  field.  The 
index  N jEg-pi  is  then  tested  and  found  to  be  positive.  The  range  of  the  new 
BH  curve  is  computed, 


HHEI  '  HDAI  (HHC,  J-l,  I  "  HHCJI** 


(69) 


Two  of  the  eleven  stored  polynomials  are  selected  for  interpolation.  This  is 
accomplished  by  taking  the  difference,  -q^IC’  between  the  beginning  value, 
HpA  j r  j  of  each  of  the  stored  polynomials  and  the  D-A -indexed  value  of 
the  flux  density,  corresponding  to  the  beginning  point  of  the  new  BH 

curve  which  is  to  be  interpolated. 


HGAIC  ‘  HFA,  IC,  1  "  HDAI  hbcji- 


(70) 


The  index  1^  runs  from  one  to  ten. 

If  Hp  is  found  to  be  positive  with  1^  equal  to  one,  it  means  that  the  beginning 
point  lies  outside  the  end  of  the  major  loop.  (This  may  be  due  to  roundoff 
error,  or  a  poor  choice  of  the  major  loop  in  the  first  place. )  In  this  case,  1^, 
is  verified  as  being  one,  and  the  major  loop  is  used  for  computing  the  flux 
density.  The  argument  of  the  polynomial  is  the  D-A -indexed  difference  between 
the  field  and  the  beginning  point  coordinate, 


HHDI  "  HDAI  (H1I  '  HHCJ1)# 


(71) 


The  index  is  set  eQua^  t0  tw°.  indicating  that  the  upper  interpolation 

curve  is  the  one  above  the  major  loop,  and  the  lower  one  is  the  major  loop. 
The  interpolation  coefficient  is  set  equal  to  zero,  so  that,  in  effect,  the  major 
loop  is  used.  The  value  of  the  polynomial  is  Hv,  and  the  flux  density  is  found 
from  Equation  (38).  The  new  value  of  the  previous  field  for  the  next  time 
increment  is  set  equal  to  the  old  present  value  by  using  Equation  (39). 


If  is  found  to  be  positive  for  a  certain  value  of  1^  between  two  and  ten, 

inclusive,  the  stored  curve  corresponding  to  this  value  is  the  upper  interpola¬ 
tion  curve.  The  curve  corresponding  to  1^-1  is  the  lower  interpolation  curve. 


If  the  difference  given  by  Equation  (70)  is  found  to  oe  negative  with  1^  equal 
to  ten,  the  index  1^,  is  set  equal  to  eleven. 


After  1^  has  been  determined  to  be  between  two  and  eleven,  inclusive,  the 
beginning  point  value  of  the  interpolation  coefficient  is  found  as  the  ratio  of 
the  appropriate  differences  in  flux  density  coordinates, 


HKAI 


-H 


GA.IC-1 


HGAIC  ’  HGA,IC-1 


(72) 


The  value  of  the  flux  density,  H^.  from  the  lower  interpolation  curve,  I^-l, 
is  found  from  its  polynomial,  using  the  range  found  from  Equation  (69)  as  the 
argument.  Similarly,  the  value,  H^,  from  the  upper  interpolation  curve, 
1^,  is  found. 


The  end  point  value  of  the  interpolation  coefficient  is  found  as  the  ratio  of  the 
appropriate  flux  density  differences, 


HKCI 


hdai  HBC,  J-l.I 
^IUI  '  HILI 


(73) 


The  index 

rCURI  =  IC’ 


(74) 


is  used  henceforth  to  indicate  the  upper  interpolation  curve. 


The  argument  of  the  polynomial  is  computed  from  Equation  (71).  This  is  used 
in  the  polynomials  representing  the  interpolation  curves.  The  value  from  the 
lower  interpolation  curve  is  designated  and  that  from  the  upper  one  is 


H 


FUI* 
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The  parabolically  sliding  interpolation  coefficient  is 


HKLI  "  HKAI  +  (HKCI  "  HKAI) 


(75) 


The  flux  density  is  computed  by  interpolation  between  the  interpolation  curves, 

HBLI  =  HDAI  [hkli  HFUI  +  ^1'hkli)  hfli]  *  '76) 

The  new  value  of  the  previous  field  for  the  next  time  interval  is  set  up  by  using 
Equation  (39). 


The  progi'am  then  returns  to  the  output  subroutine,  and  the  required  variables 
are  stored  on  tape  for  later  printout. 


Procedure  With  No  Reversal  of  H  -  With  no  reversal  of  H,  the  D-A  index  re¬ 
mains  unchanged.  Nevertheless,  the  value  of  the  field  must  be  tested  to  de¬ 
termine  whether  it  exceeds  the  end  point  of  the  present  BH  curve.  This  is 
done  by  means  of  Equations  (62)  and  (63).  The  same  procedure  is  followed 
as  in  the  previous  section.  When  inequality  (68)  is  satisfied,  a  curve  with  a 
satisfactory  end  point  has  been  found.  If,  during  this  procedure,  inequality  (64) 
was  found  to  hold  one  or  more  times  the  index  will  have  been  set 

equal  to  one.  Under  this  condition,  the  procedure  of  the  last  section  will  be 
followed,  from  Equation  (69)  onward. 


If  inequality  (68)  was  satisfied  the  first  time  that  was  tested,  the  value 

of  H  did  not  exceed  the  end  point  of  the  present  BH  curve,  so  no  new  curve  ne^d 
be  derived.  The  index  N-pggj  is  tested  and  found  equal  to  zero.  Then  the 
procedure  of  the  last  section  will  be  followed,  using  Equation  (71)  to  find  the 
argument  of  the  polynomial,  the  interpolation  curves  to  find  Hppj  and  Hppj, 
and  continuing  with  Equations  (75),  (76),  and  (39). 


4.  3  CURVE-FITTING  AND  INTERPOLATION  TECHNIQUES 

The  fitting  of  analytical  expressions  to  the  aclual  BH  curves  for  the  material,  and  the 
generation  of  an  interpolation  scheme,  proved  to  be  more  difficult  than  the  generation  of 
the  logical  scheme.  This  difficulty  was  due  to  the  shape  of  the  BH  curves,  to  the  limited 
storage  in  the  core  memory  of  the  computer,  and  to  the  necessity  of  obtaining  economical 
running  speeds.  The  last  two  factors  required  that  the  logic  be  simplified  as  much  as 
possible,  that  a  BH  curve  be  derived  from  the  least  possible  amount  of  information,  and 
that  the  minimum  amount  of  information  about  previous  BH  curves  should  be  stored. 


With  all  of  these  factors,  as  well  as  accuracy  of  the  fit,  taken  into  account,  the  final 
choice  for  the  analytical  function  was  a  seventh  degree  polynomial.  The  manner  of 
derivation  of  the  coefficients,  the  interpolation  scheme  used  in  the  digital  computer 
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program,  and  related  matters,  are  discussed  in  Section  4.  3. 1  below.  Other  methods 
which  evidently  have  merit  are  discussed  in  Section  4.  3.  2.  The  methods  which  were 
investigated,  but  found  unsatisfactory,  are  described  in  Section  4.  3.  3.  These  are  in¬ 
cluded  for  the  purpose  of  documentation  and  to  prevent  future  effort  from  being  expended 
along  these  lines  without  a  realization  of  the  shortcomings  already  found. 

4.3.1  Seventh  Degree  Polynomials 

BH  curves  for  the  rod  material  were  obtained  from  the  Allegheny  Ludlum  Steel  Corpora¬ 
tion,  and  processed  as  described  in  Section  4.  3. 1.  1  below.  The  fitting  of  the  polynomials 
and  the  interpolation  of  additional  reference  curves  are  described  in  Section  4.  3. 1.  2. 

The  interpolation  scheme  used  in  the  digital  program  is  described  in  Section  4.  3. 1.  3. 


4.  3. 1. 1  Allegheny  Ludlum  Curves 

Several  sets  of  curves  were  obtained  from  the  Allegheny  Ludlum  Steel  Corporation.  The 
set  which  was  used  is  reproduced  as  Figure  4-3.  The  curves  were  taken  for  a  ring 
sample  of  the  material,  and  so  they  must  be  corrected  for  the  demagnetization  effect 
when  applied  to  rods.  For  the  rod  dimensions  used,  the  demagnetization  factor,  D^, 
is  1.  5  x  IQ’5  oersteds  per  gauss.  The  correction  is  applied  to  the  field, 


H  =  H-  +  D..B, 
e  l  M  ’ 


(1) 


where  B  is  the  flux  density  at  the  center  of  the  red,  1L  is  the  internal  field  at  the  center 
of  the  rod,  and  He  is  the  applied  field,  which  is  external  to  the  rod.  1L  is  then  the  value 
of  the  field  which  is  plotted  as  the  abscissa  of  the  original  set  of  curves  in  Figure  4-3. 

For  the  present  purpose,  Hg  is  the  component  along  the  rod  axis  of  the  earth’s  magnetic 
field  at  the  satellite  location.  Hg  is  the  abscissa  of  the  derived  set  of  curves  in  Figure 
4-4.  The  derived  set  was  obtained  by  replotting  with  a  new  abscissa  value,  in  accordance 
with  Equation  (1).  The  values  of  He  were  actually  obtained  from  the  first  part  of  a 
digital  computer  program  which  is  described  in  Appendix  II. 

4.  3.1.2  Curve-Fitting 

The  analytical  expression  used  for  fitting  the  BH  curves  is  a  seventh-degree  polynomial. 
The  coefficients  were  derived  by  means  of  the  least  squares  method.  The  digital  com¬ 
puter  program  used  is  described  in  detail  in  Appendix 

The  argument  of  the  polynomial  is  the  translated  value  of  the  field,  that  is  the  difference 
between  the  value  of  He  at  any  point  and  the  value  at  the  beginning  point  of  the  curve 
being  fitted.  Likewise,  the  ordinates  were  translated  by  subtracting  the  beginning  point 
value  from  each  of  the  actual  values.  For  the  translated  curve,  the  beginning  point  then 
has  the  coordinates  (0,0).  The  polynomial  was  forced  to  pass  through  this  point  by  omitting 
the  constant  term  from  the  general  expression.  To  obtain  the  reference  curves  used  in 
the  digital  program,  the  flux  density  was  retranslated  by  adding  the  beginning  point  value, 
which  thus  becomes  the  constant  in  the  polynomial.  The  argument  of  the  polynomial, 
however,  remains  the  translated  value  of  the  field.  This  is  most  useful  in  the  main 
digital  program,  particularly  for  interpolation,  as  will  be  apparent  later. 
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This  procedure  was  applied  to  the  ascending  branch  of  the  major  hysteresis  loop  and  to 
the  four  second-class  curves  of  Figure  4-3.  The  second-class  curves  are  the  ascending 
ones  which  begin  at  various  points  along  the  descending  branch  of  the  major  loop.  The 
polynomials  thus  derived  are  plotted  in  Figure  4-5,  and  designated  as  curves  number 
1,  3,  5,  7,  and  9.  These  designations  correspond  to  the  reference  curve  numbers  used 
in  the  main  program.  Reference  curve  11  is  simply  a  horizontal  line,  and  is  used  only 
when  necessary  to  interpolate  a  curve  above  the  uppermost  BH  reference  curve  available. 
(Attempts  to  use  the  descending  branch  of  the  major  loop  for  this  purpose  yielded  un¬ 
satisfactory  results. ) 

Reference  curves  2,  4,  6,  8,  and  10,  also  shown  in  Figure  4-5,  were  obtained  by  inter¬ 
polating  between  the  odd-numbered  curves. 

The  interpolation  scheme  described  in  the  next  section  requires  extrapolating  each  of  the 
curves  numbered  2  through  10  to  reach  the  abscissa  of  the  end  point  of  tne  unextrapolated 
curve  immediately  below.  These  extrapolated  sections  are  shown  dashed  in  Figure  4-5. 
Each  polynomial  was  .cted  ov<  the  complete  range,  including  the  extrapolated  portion. 

The  polynomials  are  of  the  form 

8  •  i 

B  =  £  A.  HJ_1  ,  (2) 

J=1  J 

where  H  is  the  translated  field. 


It  became  apparent  in  the  course  of  the  analysis  that  the  maximum  value  of  tne  earth's 
magnetic  field  at  the  perigee  altitude  would  be  about  0.  55  oersteds.  This  value  is  greater 
than  the  0.  405  oersteds  corresponding  to  the  tip  of  the  major  loop  in  Figure  4-4.  To 
overcome  this  deficiency,  additional  curves  were  requested  from  Allegheny  Ludlum. 

These  have  not  yet  been  received.  Therefore,  the  curves  of  Figures  4-4  and  4-5  were 
adapted  to  the  requirements  of  the  program.  This  was  accomplished  by  applying  different 
stretching  factors  to  the  field  and  flux  density  scales.  The  factors  used  were  derived 
by  extrapolating  the  magnetization  curve.  This  curve  was  assumed  to  go  through  the 
tips  of  the  major  loops  of  the  2,  4,  6,  and  8  kilogauss  sets  of  Allegheny  Ludlum  curves. 
The  coordinates  of  these  points  are  tabulated  in  Table  4-1. 


TABLE  4-1.  POINTS  ON  MAGNETIZATION  CURVE 


External  Field,  H 
(Oersteds) 

0.  084 
0.  140 
0.  202 
0.405 


Flux  Density,  B 
(Kilogauss) 

2.0 
4.0 
6.  0 
8.0 


A  parabola, 

He  -  0.  436  -  0. 144  B  +  0.  0175  B2, 


(3) 
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FLUX  DENSITY  ~  K1LOGAUSS 


was  fitted  to  the  last  three  points  in  Table  4-1.  Figure  4-6  is  a  plot  of  the  magnetization 

curve  plotted  from  Table  4-1.  A  portion  of  the  extrapolating  parabola  is  shown  as  a 

dashed  curve  in  this  figure.  It  will  be  noted  that  the  abscissa  is  the  external  field,  H  . 

’  e 

It  is  easy  to  show  algebraically,  with  the  use  of  Equation  (1),  that  if  He  is  a  second-degree 
polynomial  in  B,  then  is  also.  Thus  equivalent  results  would  have  been  obtained  if  the 
extrapolation  had  been  done  in  terms  of  IT. 

To  allow  some  margin,  the  stretching  factor  for  flux  density  cas  taken  as  1.  15.  This 
brings  the  8-kilogauss  peak  value  to  9.  2  kilogauss.  The  corresponding  value  of  He 
from  Equation  (3)  is  0.  5924  oersteds.  The  ratio  of  the  former  0.  405  oersted  value  to 
this  number  is 

C  -  ^^§4  -  0.68366  .  (4) 


This  is  the  reciprocal  of  the  stretching  factor  for  the  magnetic  field  values. 


The  new  polynomial  corresponding  to  Equation  (2)  is  then 


8 


f 


(5) 


where  the  subscript  S  indicates  the  stretched  curves.  This  subscript  is  dropped  later. 
The  new  coefficients  for  the  polynomial  were  found  from 

A  Q  =  1.  15  C^"1  A.  .  (6) 

P  J 

Equation  (6)  was  applied  to  the  coefficients  of  each  of  the  eleven  reference  curves  of 
Figure  4-5.  The  results  were  used  as  the  input  reference  curves  for  the  computer  pro¬ 
gram. 

4.  3. 1.  3  Interpolation  Scheme 

Ar.  examination  of  the  curves  of  Figure  4-5  shows  that  interpolation  with  a  constant  in¬ 
terpolation  coefficient  is  not  realistic.  Curve  5  may  be  used  as  an  example.  Suppose 
that  this  curve  were  not  available,  and  it  was  required  to  obtain  it  by  interpolation.  The 
interpolation  coefficient  is  the  ratio  of  differences  between  appropriate  ordinates, 

B5-B3 

k3  =  B?-B3  ’ 

where  the  subscripts  refer  to  the  curve  numbers.  Having  the  curve  5  available,  one 
may  obtain  the  interpolation  coefficient,  k^,  as  a  function  of  the  abscissa  (translated 
magnetic  field).  This  function  is  plotted  in  Figure  4-7.  The  beginning  point  value  is 
one-half,  and  for  low  values  of  the  abscissa,  kg  remains  nearly  constant  at  this  value. 

It  then  gradually  rises  in  value.  Over  the  right-hand  portion  of  the  curve,  the  points 
seem  erratic.  This  has  been  attributed  to  instrument  errors  in  the  original  data,  in¬ 
accuracies  in  the  reading  of  the  graphs,  and  the  high  sensitivity  of  the  coefficient  to 
these  errors  in  this  region  of  the  curves. 
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The  scheme  to  be  used  for  interpolation  in  the  main  program  should  take  into  account 
the  general  shape  of  the  curve  in  Figure  4-7,  and  yet  the  scheme  must  be  fairly  simple. 
The  requirements  of  the  main  program  are  that  a  BH  curve  be  derivable  from  the  co¬ 
ordinates  of  ito  end  points,  together  with  the  use  of  the  eleven  stored  reference  curves. 

It  is  evident  that  the  interpolation  coefficient  must  generally  have  a  value  at  the  end  point 
of  the  interpolated  curve  different  from  the  value  at  the  beginning  point.  In  the  program, 
the  beginning  point  value,  H^j.  and  the  end  point  value,  Hj^j,  are  found  from  the 
equivalent  of  Equation  (7),  using  the  flux  density  at  the  beginning  or  end  point,  as  appro¬ 
priate,  in  place  of  Bg.  The  coefficient  at  ^ny  intermediate  point  is  then  taken  as  a 
parabolic  function  of  the  abscissa. 


H 


KLI  “  HKAI 


+  X2  (H 


KCI  '  HKAI)* 


(8) 


where  X  is  the  ratio  of  the  abscissa  to  its  total  range.  Therefore  X  varies  from  zero 
to  one  over  the  range  of  the  derived  curve.  It  is  easily  verified  that  Equation  (8)  satisfies 
the  beginning  and  end  point  requirements,  and  has  a  zero  slope  at  the  beginning  point. 

Such  a  function  was  considered  to  best  meet  the  requirements  of  simplicity  with  a  shape 
similar  to  that  of  Figure  4-7. 

The  scheme  was  tested  by  applying  it  to  one  of  the  Allegheny  Ludium  curves.  The  curve 
tested  is  the  longest  descending  branch  of  a  minor  loop  illustrated  in  Figure  4-4.  This 
curve  was  selected  because  it  has  a  greater  range  than  any  of  the  other  corresponding 
curves  in  its  set.  The  actual  curve  was  inverted  about  the  origin  to  make  it  ascending, 
translated,  and  replotted  in  Figure  4-8.  The  curve  obtained  by  interpolation  between 
reference  curves  4  and  5  is  shown  dashed  in  Figure  4-8  for  comparison.  The  discrepancy 
is  believed  to  be  comparable  to  the  instrument  and  other  errors. 

4.  3.  2  Alternate  Methods  Which  Appear  Applicable 

Several  methods  of  fitting  BH  curves  with  analytical  functions  and  of  interpolating  between 
such  curves  were  investigated.  One  of  the  curve-fitting  methods,  described  later,  was 
tried  and  found  successful  but  not  used,  because  the  seventh  degree  polynomial  is  simpler. 

The  problems  fall  into  two  distinct  areas,  curve-fitting  and  interpolation.  Interpolation 
methods  are  divided  into  three  classes.  The  first  of  these  is  interpolation  between  sets 
of  curves  corresponding  to  various  peak  values  (  end  points  of  the  major  hysteresis 
loop,  such  as  the  2,  4,  6,  ar.d  8  kilogauss  sets  of  curves  obtained  from  Allegheny  Ludium. 
The  second  is  interpolation  between  the  curves  in  any  such  set,  that  is,  the  generation 
of  additional  curves.  The  third  is  the  type  of  interpolation  that  is  suitable  for  use  in  the 
digital  simulation,  that  is,  a  scheme  which  will  generate  a  curve  from  the  coordinates 
of  its  end  points,  under  the  conditions  of  an  arbitrary  history,  and  using  a  limited  amount 
of  input  information  derived  from  the  original  BH  curves.  Only  the  first  two  types  of 
interpolation  are  discussed  in  this  section.  The  third  type  was  discussed  in  Section  4.  3. 1. 
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EXTERNAL  FIELD  OERSTEDS 


Figure  4-6.  Magnetization  Curve  and  Extrapolation  Parabola 


Figure  4-7.  Interpolation  Coefficient  Figure  4-8.  Comparison  of  Actual 

and  Interpolated  Curves 
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This  section  will  discuss  curve-fitting  in  general,  then  fitting  magnetic  hysteresis  curves 
in  particular,  and  finally  the  various  types  of  interpolation. 

4.  3.  2. 1  Curve-Fitting  By  Interpolation  Techniques 
The  most  familiar  form  for  fitting  a  curve 


y  -  f(x),  (l) 

is  by  means  of  polynomials  P(x).  These  are  generally  chosen  as  "interpolation”  poly¬ 
nomials,  that  is  so  as  to  pass  through  n  points  of  the  curve  (1),  generally  equally  spaced 
over  the  interval  over  which  the  fit  is  desired,  the  degree  of  the  polynomial  being  one 
less  than  the  number  of  fitted  points,  so  that  it  has  as  many  constants  (coefficients)  as 
the  number  of  points.  If  the  fit  is  not  too  good  between  the  points,  then  one  increases 
the  number  of  fitted  points  as  well  as  the  degrees  of  the  polynomial.  Occasionally  it 
happens  that  even  this  does  not  help,  and  as  the  number  of  fitted  points  increases,  the 
discrepancies  between  f(s)  and  P(x)  between  the  x.  also  increase.  This  is  likely  to 
happen  when  the  curve  (1)  has  both  flattish  and  steep  portions. 


Some  of  these  difficulties  may  be  avoided  by  using  other  forms  for  the  curve-fitting  com¬ 
ponents.  In  particular,  we  consider  approximations  of  the  form 


f(x)  A./ 
i 


£  At  RU-Xj), 


(2) 


where  x^  are  equally  spaced,  a  distance  Ax  apart  over  the  interval  in  question,  and  B 
is  a  predetermined  constant.  The  Ai  are  chosen  so  that  the  right-hand  member  of  (2) 
agrees  with  f(x)  at  the  points  x^ 

If  the  constant  B  is  small  compared  with  Ax,  then 

R(0)  =  (3) 

B 

is  large,  and  R(x)  falls  off  rapidly,  being  equal  to  l/|^B2  +  (Ax)2j  at  x  =  Ax.  Therefore 
each  coefficient  Ai  in  (2)  is  approximately  proportional  to  the  local  ordinate  f(x^).  While 
this  gives  one  an  immediate  estimate  for  the  constants  A.  and  their  variation  with  the 
shape  of  f(x),  it  has  the  disadvantage  that  between  the  points  x.  the  right  hand  member 
of  (2)  decreases  to  numerically  small  values.  Thus  the  right-hand  member  of  (2)  has 
the  appearance  shown  schematically  in  Figure  4-9.  On  the  other  hand,  if  B  is  much 
larger  than  Ax,  then  the  shape  of  R(x)  is  rather  iike  a  parabola  with  a  small  curvature, 
at  least  for  some  distance  to  each  side  of  its  maximums.  This  follows  from  the  expansion 

R(x)  =  l/B2  -  x2/B4  +  ....,  |x|<  B.  (4) 

While  there  is  no  difficulty  in  solving  the  linear  equations  for  A^,  the  variation  of  the 
coefficients  with  the  stnpe  of  f(x)  is  not  immediately  predictable. 
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Figure  4-9, 


From  the  above  it  would  appear  that  any  special  advantage  that  the  form  (2)  might  pos¬ 
sess  would  be  attained  by  having  B  of  the  same  order  of  magnitude  as  Ax.  It  is  believed 
that  rather  than  choosing  a  large  B/Ax,  it  is  preferable  to  double  the  number  of  points 
x^.  Of  course,  the  final  decision  on  whether  the  form  (2)  is  better  than  some  other 
form  of  approximating  functions  lies  in  comparing  the  discrepancies  between  the  true 
shape  and  the  approximating  shape,  at  points  between  the  fitted  points.  A  predictable 
and  systematic  variation  of  the  coefficients  A.  with  the  shape  of  f(x)  is  also  desirable, 
if  possible. 

One  property  of  Equation  (2)  (and  this  is  true  for  any  shape  R(x)  provided  that  x^  are 
equally  spaced  points),  is  that  the  matrix  of  the  coefficients  of  the  equations  to  be  solved 
for  Ai  has  the  same  elements  in  each  diagonal  parallel  to  the  main  diagonals.  The 
matrix  is  also  symmetric  about  the  main  diagonal,  if  R(x)  is  even  in  x.  The  same  state¬ 
ments  apply  to  the  inverse  matrix. 

As  an  alternative  one  may  use  the  form  (2),  but  choose  the  coefficients  in  some  other 
way  than  by  imposing  an  exact  fit  at  x^  Thus,  as  discussed  later  in  Section  4.  3.  2.  2, 
one  may  determine  the  coefficients  in  (2)  by  the  least  squares  criterion,  that  is,  so 
that  the  sum  of  the  squares  of  the  discrepancies  between  the  left  and  the  right  sides  of 
(2)  at  the  points  xi  is  least.  \s  will  be  shown,  this  gives  rise  to  linear  equations  for 
the  coefficients  A  . 

The  form  was  used  to  fit  each  of  the  curves  in  the  8-kilogauss  set  furnished  by  Allegheny 
Ludlum.  An  eleven-point  fit  was  made  by  an  IBM  7094  digital  computer  program  coded 
in  Fortran  II.  The  equations  programmed  are  listed  in  Appendix  IE.  The  runs  were 
made  for  values  of  the  ratio  of  B/Ax  equal  to  0.  25,  0.  5,  1,  2,  and  3.  Twenty-one 
points  were  taken  from  each  curve  to  be  fitted,  the  two  end  points  and  every  0.  05  of  the 
range  of  H.  The  multiples  of  0. 1  of  the  range  of  H,  together  with  the  two  end  points, 
were  the  eleven  points  used  for  curve-fitting.  The  ten  intermediate  points  were  used 
for  checking  the  goodness  of  fit.  It  was  felt  that  the  discrepancies  between  the  fitted 
expression  and  the  original  data  would  be  greatest  at  these  points,  which  lie  midway 
between  the  fitted  points.  The  errors  were  least  for  ratios  of  B/Ax  equal  to  1.  5,  2,  or 
3.  For  any  fitted  curve,  the  ratio  selected  should  be  that  which  yields  the  least  value 
of  the  maximum  error.  After  this  selection  was  made  for  each  of  the  fitted  curves, 
the  worst  error  was  380  gauss. 


A  great  deal  of  the  above  readily  carries  over  to  other  shapes  R(x)  which  have  a  maximum 
at  x  =  0;  are  symmetric  in  x,  and  decrease  to  zero  to  either  side.  Several  such  possible 
shapes  are 


l/(x‘ 


2  1 
B2) 


2/  2 
,-x  /a 


and  sech  n(x/a). 


(5) 


If  R(x)  is  chosen  as  indicated  in  Figure  4-10,  then  (2)  reduces  to  a  polygonal  approxima¬ 
tion  to  f(x)  obtained  by  drawing  straight  line  segments  between  the  ooints  of  (i)  corres¬ 
ponding  to  equally  spaced  ordinates. 

R  (X) 


Figure  4-10. 

This  requires  merely  storage  of  the  ordinates  f(xi)  and  linear  interpolation  between  them. 
It  is  too  rough  for  high  accuracy,  but  quadratic  interpolation  based  on  interpolation  using 
several  adjacent  ordinates  may  be  very  effective. 

For  1  W//M  1 

x.  =  i(Ax)  <  x  <(i+l)  Ax  =  x.  .  1  I  I  1  (6) 

i  1*1  X,  x1+1 

where  neither  x.  nor  xi+j  is  an  end  point.  The  interpolation  is  given  by  the  interpolation 
polynomial  fitting  the  values  of  f(x^)  at  x  =  x._p  x.,  x.+1  xj+2: 


where 


f(x)  =  mQ  ffxj  j)  +  ml  f(x.)  +  m2 

mQ,  mp  m2,  m3  are  given  by 

m  (x-xi>  (x-xitl)  (x-xit2) 

°  (Ax)3  1-  2-3 


f(xi*L>  *  m3 


mQ  (x)  =  - 


) 


where 


(x-Xj.,)  (x-xui)  (x-xit2) 
m  -I  -  jj 

1  (Axr  1-1-2 

m2  =  mx  (xt  +  x.+1  -  x) 

m3  =  mQ  (xt  +  x.+1  -  x) 


mt  (x) 


(9,l)(e-l)(9-2) 
1 !  2! 


(7) 


(8) 

(9) 
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If  Xj  or  j  is  an  end  point,  then  one-sided  interpolation  may  be  used  in  place  of  (7). 
Thus  inxQ<  x  <  xQ  +  Ax  Equation  (7)  holds,  but  with  one  sided  values.  Thus  for 
xt  <  x  <  x2 


f(x)  =  mQ  f(x2)  +  . . . 


m 

o 


(x-x1)(x-x2){x-x3) 
3!  (dx)3 


(10) 


These  can  also  be  expressed  in  terms  of  first,  second,  and  third  finite  differences. 

4. 3. 2. 2  Curve-Fitting  by  Least  Squares  Techniques 

Let  f(x)  be  a  given  function  of  the  independent  variable  x  over  a  proper  interval  (a,b). 
We  are  concerned  with  approximating  to  f(x)  by  means  of  a  linear  combination  of  given 
functions : 


f(x)  ~  Cjf^x)  +  C2f2(x)  +  . . .  +  Ckfk(x), 


(1) 


where  f^,  f2,  . .  are  functions  of  x  only,  and  the  Cj  are  constants  to  be  chosen  so  that 
the  "errors"  at  n  prescribed  points 


xl»  x2'  •  • »  xn 


(2) 


are  least  in  the  sense  that  the  sum  of  their  squares  is  smaller  than  for  any  other  set 
of  coefficients  Cj  in  (1). 

It  is  evident  that  unless  the  number  n  of  points  in  (2)  is  greater  than  the  number  k  of 
terms  on  the  right  of  (1),  the  errors  can  be  made  to  vanish  (provided  that  the  determinant 
J  fj(Xj)  J  does  not  vanish),  and  the  least  square  fit  reduces  to  an  interpolation  fit 

By  a  linear  change  of  variable  the  interval  (a,b)  can  be  transformed  into  the  interval  (0, 1), 
and  this  we  assume  to  be  the  case  in  the  following.  The  points  (2)  will  then  lie  in  (0, 1). 
They  may  include  one  or  both  of  the  end  points. 


The  quantity  to  be  minimized  is  thus 


Q  = 


n 

£ 

i=l 


£  W 


(3) 


Evidently,  Q  is,  in  general,  a  definite  (that  is  a  positive)  quadratic  form  in  the  coefficients 
Cj.  It  takes  on  its  minimum  at  the  values  of  Cj  which  satisfy  the  k  equations 


1  Afi. 

2  ^ 


i=l  1=1 


(4) 


Rearranging  the  terms  results  in  k  linear  equations, 

|ici  [  J?!  W  •  w]  ,<xi),i(xt);t=l---k 


(5) 
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In  these  equations  Cj  are  the  unknowns,  while  the  coefficients  of  these  unknowns  are 
sums  of  products  of  the  functions  fj  at  the  points  (2),  and  the  right-hand  constants  are 
sums  of  products  of  the  values  of  the  left-hand  member  of  (1)  by  the  values  of  the  func¬ 
tions  fj,  all  evaluated  at  respective  values  of  the  points  (2)  chosen  properly  in  each 
factor.  The  solution  of  these  equations,  when  substituted  in  the  right-hand  side  of  (1), 
yields  the  least  square  fit  for  f(x). 


An  example  is  that  in  w^ich  the  functions  f ^  are  polynomials: 

f,(x)  =  xi'1,  (6) 


Then  Equations  (5)  reduce  to 


k 

L  Ci 

J=1  J 


n 

E 

i=l 


J-1  x.1-1 

*i  xi 


E  f(x.)  x.l'l\l=l,  k 
i-1  1  1 


(7) 


The  least  square  approximation  in  general  will  not  yield  a  perfect  fit  at  the  points  (2). 

In  particular,  if  the  end  points  are  included  among  the  points  (2),  it  is  not  to  be  expected 
that  the  least  square  fit  will  yield  the  exact  values  of  f(x)  «.t  these  points  either.  If  it 
is  desired  to  obtain  an  exact  fit  at  these  points,  say  for  polynomial  fits,  then  one  proceeds 
as  follows. 


We  choose  for  the  polynomials  the  functions 

f j(x)  =  x(l-x)  Xs*1  (8) 

thus  replacing  (1)  by 

f(x)  ~  x(l-x)  [c^  +  Cj  x  +  C2  x2  +  . . .  +  Cfexk  j  =  c.  (x*+1  -  xJ+2  ^  (9) 


where,  with  a  slight  change  of  notation,  the  subscript  has  been  allowed  to  range  over  0 
as  well.  The  form  (9)  automatically  vanishes  at  the  end  points  x  ■  0  and  x  =  1. 


As  an  alternative  method,  one  may,  before  applying  least  squares,  subtract  from  f(x) 
a  linear  function  of  x  which  takes  on  at  the  end  points  the  same  values  as  f(x).  Then  one 
proceeds  to  apply  least  squares  to  the  resulting  difference.  We  shall  assume  in  the 
following  that  this  has  been  done,  and  will  use  f(x)  for  the  difference  function.  With  this 
in  mind,  Equations  (5)  apply  and,  more  explicitly,  yield 


k 

L 

j=l 


ci 


n 

L 

i=l 


x-2  (l-xt)2  xt 


n 

=  L 

i=l 


f(xt)  x.(l-x.)  x.4 


9 


-^=1,  k. 


(10) 


Another  alternative  is  to  divide  both  sides  of  (9)  by  x(l-x)  and  choose  a  polynomial  repre 
sented  by  the  bracketed  factor  in  (9)  so  that  it  fits  f(x)/x(l-x)  in  the  least  squares  sense. 
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The  resulting  coefficients  C ^  will  not  be  the  same  as  above.  There  will  also  be  the 
difficulty  of  obtaining  limits  of  indeterminate  forms  when  evaluating  the  function 


at  the  end  points. 

Least  square  methods  can  be  readily  extended  to  the  case  where  the  squares  of  the 
deviations  are  weighted  with  weights  w.  which  are  prescribed  with  the  position  x4. 

It  is  believed  that,  as  compared  with  an  interpolation  polynomial  which  takes  on  the  values 
f(xj)  at  x-,  the  least  square  polynomials  have  an  over-all  smaller  error  over  the  whole 
interval,  say  when  x.  are  equally  spaced.  The  error  of  the  interpolation  polynomial  can 
be  expressed  in  terms  of  the  (n+l)th  derivative  of  f(x);  no  estimate  of  the  error  of  the 
least  square  approximation  polynomial  is  known. 

4. 3. 2. 3  Curve-Fitting  of  BH  Curves  by  Inversion 

Examination  of  the  hysteresis  curves  and  of  the  curves  obtained  by  reversing  the  sign  of 
dH/dt,  at  Hr  -  HR  . ,  HR  •  •  •  shows  that,  after  reversal,  the  new  curve  resembles 
the  ’’inverted  image"  of  the  preceding  curve  segment,  obtained  by  turning  the  latter 
through  180°,  about  the  midpoint  of  the  line  segment  joining  its  end  points.  This  is 
exactly  true  for  the  ascending  branch  of  the  major  loop,  which  is  obtained  from  the 
descending  one  by  rotating  the  latter  180°  about  the  origin  H=0,  B=0,  which  is  the  mid¬ 
point  of  the  line  segment  joining  Hm,  Bm  to  -Hm  'Bm*  It  is  also  true  within  very 
high  accuracy  for  any  rank  curve  and  the  curve  of  the  following  rank,  if  HR  j  and  HR  i+1 
are  close,  whereupon  the  curve  segment  between  these  reversal  points  is  a  parabola, 
and  so  is  its  inverted  image.  For  the  curves  of  rank  i,  i  >  1,  if  the  end  points  HR  ., 

Hr  are  far  apart,  and  Bm  is  large  enough  so  that  the  hysteresis  curve  has  a  point 
of  inflection  and  appreciable  reversed  curvature,  it  is  true  only  in  the  neighborhood  of 
the  reversal  point.  Thus  on  Figure  4-11  if,  after  descending  along  the  main  hysteresis 
curve  from  the  vertex  A  to  the  point  B,  H  starts  increasing,  then  the  path  followed  is 
BCA.  The  inverted  path  is  BC’  DA.  Near  B  the  two  paths  lie  close  to  each  other,  but 
they  deviate  for  larger  H. 

The  curve  BCA  passes  through  A,  the  previous  reversal  point,  and  this  is  in  accordance 
with  the  general  rule  for  the  curves  of  rank  i,  that  they  aim  for  the  previous  reversal 
point,  and  if  continued  beyond  (for  i  >  1),  follow  along  the  curve  of  rank  (i-2). 

The  inverted  curve  BC'DA  also  passes  through  the  previous  reversal  point  A,  since  in 
rotating  through  180°  about  the  midpoint  the  end  points  are  interchanged. 

If  a  curve  has  the  equation 

B  =  f(H),  HR1  <  H  <  Hm,  (1) 
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then  the  Inverted  image  between  these  points  is  given  by 

H  +  H’  =  Ha  +  HR1,  B  +  B*  =  f(H4)  +  f(HR1),  (2) 

where  half  of  the  right  hand  members  yield  the  coordinates  of  the  midpoint  of  the  straight 
line  segment  between  the  end  points,  and  H’,  B'  denote  the  point  which  is  obtained  by 
inverting  the  point  H,  B.  Hence,  in  equation  form,  the  inverted  curve  is  given  by 

B’  =  f(HA)  +  f(HR1)  -  f(H)  =  where  H  =  HA  +  HRJ  -  H\  (3) 

Eliminating  H,  and  then  dropping  the  primes, 

B  =  f{HA)  +  f(HR1)  -  f(HA  +  HR1  -H).  (4) 

With  proper  changes  of  notation  the  form  (3)  ca  be  used  to  derive  the  curves  of  rank  i 
from  those  of  rank  i-1,  where  the  latter  are  obtainable  by  inversion  of  the  preceding  ones. 

Where  an  inverted  curve  differs  appreciably  from  the  reversed  curve,  the  correction  to 
be  added  to  the  right-hand  member  of  (3)  (or  of  the  corresponding  equation  for  the  curve 
of  rank  i)  must  be  a  function  of  H  which  vanishes  at  the  two  end  points.  This  function 
increases  numerically  very  slowly  at  B  on  Figure  4-11,  or  near  its  first  reversal  point, 
while  it  approaches  its  following  reversal  point  more  steeply.  We  proceed  to  derive 
several  forms  of  functions  having  this  property.  For  definiteness,  we  change  the  inde¬ 
pendent  variable  to  the  normalized  variable  x,  and  suppose  that  the  roots  are  at  x=0, 
x=l,  with  x=0  as  the  root  with  the  slow  numerical  increase. 

A  simple  function  having  the  desired  properties  is  given  by 

y  =  C  x2  (1-x).  (5) 
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This  has  a  vanishing  slope  at  x=0,  and  would  do  only  if  tho  two  curves  (the  reversed  and 
inverted  curves)  are  iangent  to  each  other  at  the  reversal  point. 

Another  possible  function  is  given  by 

1 

y  =  C  x  (e"^x  -  e"**)  ,  £  >  0,  (6) 

where  0  and  C  are  constants.  The  ratio  of  slopes  at  x=0  and  at  x=l  is  numerically  equal 
to 

( e 0  -  1):  £  .  (7) 

I 

This  ratio  can  be  made  as  large  as  desired  by  choosing  £  large  enough.  By  replacing 
x  by  (1-x)  in  (6),  the  small  slope  can  be  made  to  occur  at  x=0,  and  the  large  one  at  x=l. 

1 

Additional  functions  with  the  desired  properties  can  be  obtained  by  starting  with  a  func¬ 
tion  like 

y  =  f(x)  =  Jn(x),  n  =  2,3,.. ,  (8) 

whose  appearance  between  x=0  and  the  first  positive  root  is  as  shown  schematically  in 
Figure  4-12,  with  a  zero  slope  at  x=0.  To  obtain  from  (8)  a  function  with  a  small  positive 
slope  at  the  left  end,  one  cuts  the  curve  v.  h  a  line 

y  =  « ,  0) 

obtaining  two  roots  x=Xj  near  zero,  x=X£  -vith  a  small  slope  at  x^.  By  replacing  x  by 

t  =  (x-x^Axg-Xj)  (10) 


one  shifts  the  roots  to  ?  =  0,  |  =  1. 


Figure  4-12. 


The  solution  of  (8),  (9)  may  have  to  be  carried  out  by  a  trial-and-error  or  by  a  successive 
improvement  method.  However,  one  may  choose  Xg  as  a  value  in  the  table  and  near  the 
root.  Then  c  *  f^)  s  known,  and  presumably  small.  The  value  Xj  can  then  be  de¬ 
termined  by  using  the  first  terra  of  the  power  series  expansion, 


Jn«  =  fn 


n 


2  n! 


(1  -...)=  t. 


(11) 
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The  precise  f(x)  to  be  used  for  AB,  for  the  family  of  2nd  rank  (H,B)  curves,  and  the 
variation  of  f(x)  with  and  with  Hm,  requires  a  great  deal  of  trial  and  error. 

4. 3. 2. 4  Differential  Equation  for  the  Second-Rank  Hysteresis  Curves 

A  family  of  curves  such  that  one  and  only  one  curve  passes  through  every  point  of  a 
certain  region  of  the  (x,  y)  plane,  can  be  described  by  means  of  its  differential  equation 

s!e  =  f(x,,),  (i) 

where  F(x,y)  is  the  slope  of  the  curve  passing  through  (x,y).  This  equation  yields  the 
curves  themselves  only  if  one  carries  out  an  integration  of  (1),  either  analytically,  or 
numerically. 

In  the  following,  an  attempt  will  be  made  to  obtain  a  differential  equation  of  the  form 

^-=F(H,B).  (2) 


for  representing  the  curves  of  rank  2,  that  is  the  curves  obtained  after  one  has  descended 
along  a  particular  hysteresis  curve  to  a  point  H^,  then  increased  H  to  Hm,  the  value  of 
H  at  the  vertex  of  the  hysteresis  curve.  The  starting  points  of  these  curves  are  the 
points  of  the  descending  branch  of  the  hysteresis  curve.  These  curves  are  shown  in 
Figure  4-13. 


Examination  of  the  curves  of  Figure  4-13  shows  that  along  any  horizontal  straight  line 

B  =  constant  (3) 

Theslopes  increase  from  a  smali  value  m^  to  the  value  m2  which  is  the  slope  of  the 
ascending  branch  of  the  hysteresis  curve  through  the  point  with  the  constant  B  in  question. 
The  value  m^  it  follows  from  the  inversion  properties  explained  in  Section  4. 3. 2. 3, 
has  the  same  slope  as  the  slope  of  the  descending  hysteresis  curve  at  its  vertex  (Hm,Bm). 


4-39 


Thus  aij  is  independent  of  the  right  hand  constant  in  (3),  depending  only  upon  Hm.  In 
fact,  on  the  Allegheny  Ludlum  (H,B)-curves  nij  appears  to  be  zero,  or  even  slightly 
negative.  This,  howe.er,  can  never  actually  be  the  case  with  any  magnetic  material, 
and  must  be  due  to  an  instrument  error. 

If,  for  simplicity,  we  do  assume  that  nij  vanishes,  then  a  possible  simple  form  for  (2) 
is  given  by 

ffi  =  6  m2  ;  <4) 

where  §  is  the  fractional  distance  along  the  line  (3)  frc*  i  the  left,  to  the  total  length  of 
this  line: 


e  =  (h  -  hd)/(ha  -  hd), 


(5) 


where 

Ha  =  Ha(B),  Hd  =  Hd(B),  (6) 

are  the  equations  of  the  ascending  and  the  descending  branches  of  the  hysteresis  curves 
respectively.  It  follows  that 

m2  =  l/dHA/dB  =  1/H^(B).  (7) 


Substitution  from  (5),  (7)  into  (4)  yields  the  differential  equation 

dB  H  '  HD<B>  1 

dH  -  Ha(B)-Hd(B)  ha'(B) 

If  mj  is  not  negligible,  then  (4)  must  be  modified,  say,  into 
dH  =  (1’6)  m!  +  6  m2 


(8) 


(9) 


Utilizing  (5)  and  (7)  one  adds  to  the  right-hand  side  of  (8) 


m. 


Ha(B)  -  H 


(10) 


By  utilizing  the  inversion  properties  of  the  two  hysteresis  branches  one  may  express  Hp 
in  terms  of  HA,  as  follows: 

Hd(B)  =  -Ha(-B)  (11) 


Whether  (8)  or  the  modified  form  with  (10)  added  on  fits  the  curves  of  Figure  4-13  can 
be  determined  by  a  numerical  comparison  of  the  actual  slopes  of  the  curves  with  those 
computed  from  the  equations,  or  better  yet,  by  integrating  these  differential  equations 
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and  comparing  the  resulting  curves  with  the  curves  of  Figure  4-13.  The  analytic 
integration  off-hand  appears  difficult.  The  numerical  integration,  if  it  is  used,  has 
to  be  carried  out  with  high  accuracy. 

Various  changes  in  (9)  are  possible  and  may  be  desirable  if  (9)  does  not  yield  a  good  fit. 
Thus,  one  may  replace  0  by  any  function  of  6,  g(0)  which  has  the  property  that  it  in¬ 
creases  with  0,  vanishes  for  0  =  0,  and  is  equal  to  1  for  0  =  1.  Moreover,  g(0)  may  vary 
with  the  height  of  the  line  of  constant  B. 

A  further  alternative  consists  in  fitting  slopes  along  vertical,  rather  than  horizontal, lines 
in  Figure  4-13.  Along  a  vertical  line,  it  will  be  seen  that,  as  B  increases,  the  slope 
starts  with  the  slope  m2  of  the  ascending  branch  of  the  hysteresis  curve  for  the  constant 
value  of  H  in  question  and  decreases  to  m^  at  the  intersection  with  the  descending  branch. 
One  may  use  (9),  but  0  defined  by 

0  =  (B-Ba)/(Bd-Ba),  (12) 


where  BA,  B^  represent  the  ascending  and  descending  branches,  with  B  expressed  in 
terms  of  H,  and 


m2  =  dBA/dH  =  Ba'(H) 


There  results 


dB  Ba’(H)[B-Ba(H)J  Bd(H)-B 

air  =  -  ba(K)  +  mi  bd(h)  -  w/ 


This  is  linear,  but  not  homogeneous,  in  B,  with  coefficients  which  are  functions  of  the 
independent  variable  H.  Such  a  differential  equation  can  be  integrated  by  means  of 
quadratures,  and  for  this  reason  (14)  would  be  preferable  to  (8)  and  its  two  consequent 
equations.  However,  the  linear  variation  of  slope  with  distance  appears  to  be  question¬ 
able  along  vertical  lines  in  Figure  4-13,  Hence  replacement  of  0  by  g(6),  a  function  of 
0 ,  may  be  necessary.  This  would  lose  the  linearity  property  of  the  differential  equation. 

While  the  use  of  differential  equations  in  describing  the  (B,  H)-curves  is  an  added  com¬ 
plication,  to  be  avoided  if  the  explicit  equations  for  these  curves  can  be  obtained,  it  may 
be  said,  that,  from  the  point  of  view  of  integrating  for  the  motion  of  a  satellite  and  its 
orientation  along  its  orbit,  the  number  of  degrees  of  freedom  is  already  high,  and  the 
use  of  (8)  or  (14)  only  increases  by  1  the  order  of  the  system.  What  is  a  more  serious 
item  is  the  fact  that  the  differential  equations  proposed  have  only  been  considered  for  the 
second  rank  (H,  B)  curves.  For  third,  and  higher  rank  curves,  which  depend  upon 
several  preceding  reversal  values  of  H  as  well,  the  form  (2)  is  insufficient  and  either  a 
higher  order  system  is  needed  or  a  way  of  describing  the  changes  in  F  upon  reversal 
(of  the  sign  of  dH/dt). 


j 
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It  may  be  added  in  conclusion  that  from  a  scientific  point  of  view,  it  would  be  very  de¬ 
sirable  if  a  way  of  describing  magnetic  properties  of  ferromagnetic  materials  based  on 
some  theory,  were  worked  out  by  physicists,  from  which,  starting  with  the  established 
tenet  that  elementary  magnets  have  a  constant  moment,  and  all  the  magnetic  properties 
are  due  to  the  orientation  of  these  elementary  (atomic  or  molecular)  magnets,  one 
could  predict  the  magnetic  behavior  of  a  magnetic  material,  once  its  chemical  and 
metallurgical  structure  is  known.  The  work  by  Ewing  of  a  model  based  on  a  collection 
of  regularly  arranged  collection  of  small  magnets,  raised  hopes  that  this  may  be  possi¬ 
ble,  but  to  date  the  task  has  proved  too  difficult,  though  the  nature  of  the  elementary 
magnets  has  indeed  been  clarified  in  terms  of  quantum  theory,  Bohr  magnetons,  and 
electron  orbital  spin.  It  is  possible,  that  when  this  theory  is  worked  out,  it  will  be  in 
the  form  of  differential  equations  resembling  (2).  A  review  of  Neel's  theory  and  other 
theoretical  attempts  in  the  booklet  by  G.  F.  Brown  shows  that  theory  is  still  a  long  way 
from  theoretical  explanation  of  properties  of  magnetic  materials.  Possibly  some  semi- 
empirical  theory,  with  the  inclusion  of  some  sort  of  "frictional"  energy  dissipation 
accompanying  (H,  B)  changes  might  be  developed  in  the  meantime. 

4.  3.  2. 5  Interpolation  of  Second-Rank  BH  Curves 

For  a  given  hysteresis  loop  between  the  points  (Hm,  Bm)  and  (-Hm,  -Bm),  second  rank 
(H,  B)  curves  are  available,  proceeding  from  five  points  of  the  descending  branch 
corresponding  to  values  of  B^  differing  by 

B  *  I  Bm  <l> 

and  corresponding  to 

B  =  Bm,  (3/5)  Bm>  (1/5)  Bm,  -(1/5)  Bm,  -(3/5)  Bm,  -  Bm>  (2) 

where  the  last  curve  is  the  ascending  branch  of  the  hysteresis  loop  and  forms  a  limiting 
case  of  the  second  rank  curves.  These  curves  start  at  va’.ues  of  H  denoted  by  HR  (cor¬ 
responding  to  reversal  of  sign  of  dH/dt)  and  all  proceed  to  the  point  (Hm,  Bm). 

The  question  is  that  of  interpolation  between  the  above  curves. 

One  method  is  that  of  linear"slanting"  interpolation.  This  consists  in  replacing  the 
descending  hysteresis  curve  by  straight  lines  between  adjacent  starting  points  of  the 
second  rank  curves,  as  shown  schematically  on  Figure  4-14  at  PQPp  This  line  segment 
PgPj  is  continued  past  Pj  until  it  meets  the  vertical  line  H  =  Hm  at  a  point  Qp  Then 
many  line  segments  are  drawn  between  the  second  rank  curves  through  Pp  P<2,  all 
proceeding  from  Qp  The  interpolated  curves  are  obtained  by  dividing  these  segments 
in  a  constant  ratio.  On  Figure  4-14  the  mid-points  on  five  such  line  segements  have 
been  indicated.  A  curve  passing  through  these  points  is  the  interpolated  curve  corres¬ 
ponding  to  an  initial  poinl  Pj  5  corresponding  to  midway  between  of  Pj  and  Pg. 
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Figure  4-14. 


It  also  corresponds  to  a  B-value  midway  between  those  of  Pj  and  Pg.  Thus  the  starting 
point  does  not  lie  on  the  descending  branch  proper,  the  discrepancy  being  large  where 
the  curvature  of  the  hysteresis  curve  is  appreciable. 


Each  segment  PqPj,  PiP2’  *  *  *  cuts  1116  vertica*  through  PQ  in  a  separate  point  Qo(-Pq), 
Qj,  Q2,  •  •  •  The  discrepancy  in  the  top-most  interval  is  most  severe. 


An  improvement  in  the  above  procedure  consists  in  replotting  the  original  second  rank 
curves  so  that  they  start  at  points  on  a  vertical  axis  and  stretching  each  one  horizontally 
in  a  proper  ratio  so  that  they  end  in  the  same  point.  The  resulting  diagram  is  shown 
schematically  in  Figure  4-15,  where  the  B  axis  has  the  same  scale  as  in  Figure  4-14, 
and  the  variable  x  is  given  by 


x 


H  -  H 


R 


H  -  Hr 
m  K 


(3) 


and  runs  from  0  to  1  for  each  second  rank  curve.  The  ratio  of  scales  is  given  by  the 
factor 

1  -  Hm  -  «R  <4> 


and  this  is  plotted  on  Figure  4-16. 

The  second  rank  curves  on  Figure  4-14  are  obtained  by  multiplying  each  curve  of  Figure 
4-15  by  the  factor  f  of  Figure  4-16  lying  at  the  same  B  as  the  starting  point  of  Figure  4-16 
and  transposing  it  to  Figure  4-14  parallel  to  the  H  axis. 

Figure  4-16  has  the  same  shape  as  the  upper  descending  branch  of  the  hysteresis  curve 
of  Figure  4-14  except  for  bein^  reflected  in  a  vertical  axis.  If  f  were  plotted  to  the  left, 
the  curve  of  Figure  4-16  would  be  congruent  to  that  of  Figure  4-14. 
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On  Figure  4-15  the  interpolation  is  curved  out  along  vertical  segments,  as  indicated 
schematically  between  the  curves  through  Pj  and  Pg.  Again,  the  midway  value  is  shown, 
as  well  as  the  curve  (broken)  passing  through  these  midway  values.  After  this  curve  is 
obtained,  it  is  stretched  horizontally  by  the  value  of  f  corresponding  to  the  initial  B. 

The  advantage  of  the  above  method  over  that  of  Figure  4-14  is  that  one  may  fit  f(B)  in 
Figure  4-16  as  a  function  of  B  by  means  of  French  curves  or  by  means  of  an  analytic 
expression,  say,  passing  through  as  many  points  as  desired,  since  this  curve,  as  indi¬ 
cated  above,  is  derivable  from  the  (descending  branch  of  the)  hysteresis  curve. 

If  the  curve  of  Figure  4-16  is  replaced  by  straight  line  segments  through  the  values  (2) 
of  B,  then  the  method  of  Figure  4-15  (and  the  modified  Figure  4-16)  becomes  identical 
with  that  of  Figure  4-14.  This  procedure  is  not  recommended. 

The  curve  fitting  of  Figure  4-16 


f  =  f(B),  Hm  -  Hr  =  f(B),  Hr  =  Hm  -  f(B),  (5) 

should  yield  Hm  as  a  function  of  B  and  not  the  other  way  around.  This  is  advisable  in 
addition  to  the  fitting  of  the  hysteresis  branch  in  the  form 

B  =  function  of  H  =  g(H)  (6) 

where  g  is  the  function  inverse  to  the  Hm  -f(B).  It  is,  of  course,  possible  to  solve  (6) 
by  a  successive  approximation  method  or  trial  and  error  for  the  H  values  corresponding 
to  prescribed  B-values,  then  identifying  H  with  HR  in  (5),  obtain  f  =  Hm  -  H. 
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The  hysteresis  curves  of  rank  one,  and  the  curves  of  rank  2  are  available  only  from 
recording  pen  tracings  of  the  magnetic  instrument,  and  exhibit  "wobbles"  and  other 
irregularities  that  have  to  be  corrected  graphically.  It  is  evident  that  the  aid  of  a  good 

draftsman  is  needed  to  draw  the  curves  of  Figure  4-14.  When  this  is  done  on  graph 
paper  with  small  squares,  many  points  of  Figure  4-14  can  be  read  off  and  the  values  of 

B  and  H  tabulated  for  the  descending  hysteresis  branch  and  the  four  second  rank  curves. 
From  the  H  values  one  can  calculate  x  values  and  plot  the  curves  of  Figure  4-15. 

It  is  planned  to  fit  the  curves  of  Figure  4-14  (first  and  second  rank  curves)  with  equations 
of  the  form 


B  =  f(H,  Hr)  or  B  =  f(H,  BR)  (7) 

From  these,  analytic  forms  can  be  obtained  for  the  curves  of  Figure  4-15,  giving  B  as 
functions  of  x  (and  the  starting  point  Bp'.  Suppose  that  the  Wm  of  the  equations  is 

B=  7  Ak(BR)fR(X)  .  0<  x<  1,  (8) 

k 

where  the  shapes  of  fk  are  given  and  are  the  same  for  all  the  five  values  (2)  and  Ak  are 
constants  depending  on  the  curve.  For  instance,  if  fk(X)  are  powers  of  x,  the  equations 
(8)  constitute  a  polynomial  fit. 

Then  the  interpolation  between  the  curves  of  Figure  4- 15  can  be  carried  out  by  applying 
it  to  the  coefficients  Ak.  Thus,  between  the  curves  through  Pj  and  P2 

B  [(6B1  +  1-S  B2),  xl 

r  i  (9) 

=  z  [flAk(B1)  +  (l-e)Ak(B2)J  fk(x),  0  <  e  <  1. 

If,  for  definition’s  sake,  there  are  10  values  of  k  in  (8)  then  a  table  of  50  coefficients 
Ak  (Bf)  is  all  that  is  required,  since  (9)  can  be  used  to  obtain  the  equation  of  any  inter¬ 
mediate  curve  in  terms  of  the  two  curves  to  each  side. 

In  principle,  it  is  possible  to  interpolate  Ak  not  linearly  as  in  (9),  but  by  using  a  Lagrange 
interpolation  value  between  three  adjacent  curves  (or  values  of  BR). 

While  we  have  assumed  that  only  second  rank  curves  and  one  hysteresis  curve  are  available, 
it  is  evident  that  the  technique  is  applicable  if  a  larger  number  of  curves,  i.  e. ,  of  higher 
rank,  were  available. 

A  check  should  be  made  between  the  interpolated  or  predicted  curve  and  a  measured 
curve,  specially  taken  for  this  purpose.  If  the  check  is  not  too  good,  the  more  complex 
interpolation  (using  say  three  adjacent  curves,  or  preferably  a  total  of  nine  secondaries 
plus  one  primary,  ana  interpolation  of  the  form  (9))  should  be  resorted  to. 
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4.  3. 2. 6  Interpolation  Between  Major  Hysteresis  Curves 

In  Section  4. 3. 2. 5  we  discussed  the  interpolation  of  the  second  rank  BH  curves  indicated 

in  Figure  4-14.  The  procedure  recommended  there  is  shown  in  Figures  4-15  and  4-16. 

It  is  evident  that  this  procedure  can  be  extended  to  other  one-parameter  families  of 

curves.  As  an  example,  we  now  consider  the  interpolation  between  the  main  hysteresis 

curves,  or  the  curves  of  rank  1.  The  latter  are  available  for  4  values  of  B  , 

m 

Bm  =  2000»  4000’  6000»  8000  S31188  (D 

We  consider  only  the  descending  branches,  indicated  schematically  in  Figure  4-17  with 
arrows.  The  upper  end  points  of  these  curves  correspond  to  equally  spaced  values  of  B, 
displayed  in  equation  (1). 

The  interpolation  may  be  carried  out  by  introducing  a  change  of  scale  for  the  horizontal 
coordinate,  H: 

x  =  H/Hm,  (2) 

and  replotting  these  curves  in  the  x,  B  plane,  where  they  will  all  proceed  from  x=l  to 
x=-l.  They  are  shown  schematically  in  Figure  4-  18L 


The  interpolation  is  carried  out  in  Figure  4-18  by  means  of  vertical  segments  joining 
adjacent  curves.  Thus,  corresponding  to  Bm  =  7500,  these  segments  are  divided  in 
the  ratio  3:1  with  the  shorter  segment  abutting  on  the  Bm  =  8000  curve.  Where  the  two 
curves  intersect,  the  interpolated  point  will  coincide  with  the  point  of  intersection,  while 
to  each  side  of  this  intersection  the  direction  of  the  segment  proceeding  from  the  6000 
curve  toward  the  8000  curve  will  have  opposite  directions. 

After  a  sufficient  number  of  vertical  segments  and  interpolated  points  have  been  obtained, 
the  complete  curve  corresponding  to  Bm  =  7500  is  drawn  through  them. 

To  transfer  the  interpolated  curve  to  the  (BH)-plane  of  Figure  4-17,  it  is  necessary  to 
stretch  it  horizontally  by  the  scale  factor 
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corresponding  to  the  Bm  value  in  question.  N^w  the  Hm  value  can  be  read  off  from 
Figure  4- 17  by  drawing  the  curve  through  the  vertices  of  the  hysteresis  loops,  shown  as 
a  broken  line  in  Figure  4- 18,  and  described  by 


B  =  B  (H  ), 
m  m  m  ’ 

or  preferably  by 

H  =  H  (B  ), 
m  m  m  ’ 


(4) 

(5) 


as  the  abscissa  corresponding  to  the  desired  Bm. 

It  is  also  possible  to  effect  two  changes  of  scale,  one  horizontal,  as  above,  the  other 
vertical,  so  that  the  height  of  the  transformed  curve  is  also  fixed.  Then  all  the  vertices 
of  the  transformed  curves  pass  through  two  fixed  end  points.  The  interpolation  can  then 
be  carried  out  either  horizontally  or  vertically.  Aside  from  saving  the  extra  change  of 
scale,  Figure  4- 18  has  the  advantage  that  the  segment  of  the  x-axis  between  x=-l  and  x=+l 
can  be  included  among  the  family  of  curves  (it  represents  the  origin  on  Figure  4-18), 
and  used  to  interpolate  for  values  of  Bm  less  than  2000. 

As  mentioned  in  Section  4. 3. 2. 5,  it  is  also  possible  to  use  more  sophisticated  interpola¬ 
tion  by  utilizing  more  than  the  two  adjacent  curves  in  Figure  4-18. 


In  equation  form,  with  rather  cumbersome  notation,  using  only  linear  interpolation,  let 
the  hysteresis  curves  be  given  by 


B  =  B(H,Bm), 


B  =  B 
m  m,  l’ 


i  =  1,  2,  3,  4. 


(6) 


Then  the  curves  of  Figure  4-18  are  given  by 


B  -  B(H  x,  B)  =f(x,Bj. 
m  ’  m  ’  m 


(7) 


The  interpolated  curve,  corresponding  to  the  value 

B  =  (1-6)  B  .  +  0  B  .  . ,  0  <  6  <  1 

m  m,i  m,i+l’ 

is  given  by 

B  =  (!-e)  f(x,Bmf  4)  +  «  f(x»Bm>i+1)* 


(8) 

(9) 


The  interpolated  curve  on  Figure  4-17  is  represented  by 

B  =  (1-e)  f(H/H  ,  B  .)  +  a  f(H/H  ,  Bm  .  .), 
m’  m,i  m’  m,i+l  ’ 


(10) 


where  Hm  is  obtained  from  Equation  (5)  as  the  value  of  Hm  corresponding  to  (8). 
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In  implementing  the  above  on  computing  machines,  suppose  that  polynomial  approxima¬ 
tions  are  used  to  fit  the  BH  curves.  Starting  with  polynomials  of  a  certain  degree  for 
the  (four)  curves  f  in  Equation  (7)  over  the  interval 

-1  <  x  <  +1  (11) 

and  assuming  a  value  of  8  in  (8)  one  carries  the  interpolation  (9)  by  interpolating  the 
respective  coefficients  of  the  polynomials.  If  the  right-hand  member  of  (5)  is  likewise 
approximated  by  means  of  a  polynomial  in  Bm,  then  its  value  for  Bm  given  by  (8)  is 
computed.  On  the  other  hand,  if  only  its  inverse  function,  represented  by  Equation  (4), 
has  been  approximated  (say  by  means  of  a  polynomial  in  Hm),  then  it  is  necessary  to 
solve  for  a  root  of  (4)  by  some  trial -and-error  or  successive  approximation  method. 

Since  H  is  known  as  a  function  of  time,  the  representation  (7)  for  the  (BH)-curves  is 
natural.  Since  this  involves  Bm,  because  the  magnetic  data  have  been  taken  for  specified 
(and  equidistant)  values  of  Bm,  equation  (4)  must  be  used  first,  assuming  that  Hm  will 
be  prescribed  with  the  satellite  motion.  Then  6  is  determined  from  (8)  and  the  procedure 
continued  as  outlined  above. 

While  the  procedure  has  been  described  as  an  interpolation  between  the  (first  rank) 
hysteresis  curves  proper,  the  method  is  obviously  applicable  to  any  one-parameter  family 
of  curves,  for  instance  to  the  2nd  rank  curves  described  in  Section  4. 3.  2. 5,  and  repre¬ 
senting  the  B  values  after  one  has  proceeded  a  certain  distance  along  a  particular  hysteresis 
descending  branch,  then  changing  the  sign  of  dH/dt  by  letting  H  increase  again  after  it 
has  reached  a  minimum  at  HR.  If  now  H  increases  continuously,  until  H  has  attained 
the  upper  vertex  of  the  previous  hysteresis  curve,  then  starts  to  decrease  again,  the 
original  (descending)  branch  of  the  hysteresis  curve  is  followed.  If,  however,  H  increase 
to  a  value  2  less  than  Hm,  then  starts  decreasing,  then  one  follows  a  3rd  rank  curve. 

The  second  rank  curves  involve  two  parameters,  in  addition  to  H: 


B  =  B(H,Bm,HR), 


and  are  available  for  a  mesh  of  Bm,  HD  values.  The  third  rank  curves  are  of  the  form 

m  rv 

B  =  B(H,Bm,HR1,HR2).  (13) 

They  evidently  involve  more  storage  space,  and  must  be  based  on  time-consuming  taking 
of  data.  After  H  has  attained  the  maximum  at  HR  2,  it  decreases,  following  (13),  to  a 
certain  minimum  HR  g,  then  starts  increasing  again.  If  HR  g  is  greater  than  HR^,  then 
a  4th  rank  curve  is  followed,  for  which  even  less  data  is  available.  On  the  other  hand, 
“HR.  g  is  less  than  Hf  j,  then  the  last  portion  of  (13)  coincides  with  the  (6),  and  further 
reversals  of  sign  of  dH/dt  can  be  followed  for  a  while  with  the  curves  (12),  (13). 

The  higher  and  higher  rank  curves,  which  are  presumably  required  if  H  oscillates  be¬ 
tween  maxima,  and  minima,  both  of  which  are  decreasing  numerically  and  approaching 
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zero,  can  be  obtained  after  a  while  by  the  inversion  method  described  in  Section  4, 3. 2. 3. 
As  the  interval  between  successive  HR  values  decreases,  these  curves  become  Rayleigh 
parabolas. 

Opinion  seems  to  be  divided  on  what  happens  if  ^  proceeds  to  a  value  larger  than  Hrj. 
In  the  absence  of  complete  data  on  this  point,  the  recommendation  is  to  proceed  along  the 
curve  (4)  beyond  the  vertex  point  of  the  original  hysteresis  curve,  and  then  proceed  along 
an  outer  hysteresis  curve  after  H  starts  decreasing  again.  This  may  be  somewhat  in 
error,  since  the  hysteresis  curves  are  generally  taken  only  after  several  reversals  of 
dH/dt  occurring  at  *  Hm.  The  "settling  down"  to  the  final  hysteresis  curve  occurs 
more  quickly  if  H  is  changing  at  a  slow  rate.  F  om  this  it  may  be  presumed  that  the 
statement  about  picking  up  the  new  hysteresis  curve  may  be  true  for  slow-swinging 
periods  (about  a  half  hour  or  more). 

4. 3. 3  Other  Curve-Fitting  Functions  Tried  and  Abandoned 

Several  methods  of  fitting  analytical  expressions  to  the  magnetic  hysteresis  curves  were 
attempted,  brought  to  various  stages  of  completion,  and  abandoned.  The  functional  forms 
tried  were  obtained  from  the  literature  or  were  modifications  of  such  functions.  Each 
of  the  forms  is  discussed,  and,  in  most  cases,  its  advantages  and  disadvantages  are 
discussed.  In  every  case,  the  reason  for  abandoning  the  particular  form  is  given.  This 
section  serves  not  only  to  document  the  work  done,  but  also  to  prevent  consideration  of 
using  any  one  of  these  forms  in  the  future,  without  due  regard  for  the  reason  that  it  was 
rejected  here. 

4.  3. 3. 1  Rayleigh  Loops 

Rayleigh  loops  are  parabolas,  for  which  B  is  a  second  degree  polynomial  in  H.  These 
are  discussed  on  pages  490-491  of  Reference  2.  The  formulas  are  given  for  loops  whose 
end  points  are  symmetrical  with  respect  to  the  origin,  but  the  loops  are  easily  translated 
to  fit  other  situations. 

The  use  of  a  parabola  allows  the  choice  of  three  coefficients.  These  are  usually  forced 
to  fit  the  two  end  points  and  the  initial  slope  of  the  curve.  This  initial  slope  may  be 
chosen  to  be  that  of  the  initial  magnetization  curve,  or  it  may  be  a  function  of  the  flux 
density  at  which  the  reversal  occurred.  Alternatively,  the  three  coefficients  may  be 
chosen  to  fit  the  two  end  points  and  the  height,  w,  of  a  symmetrical  loop,  as  illustrated 
in  Figure  4-19.  The  attempt  to  fit  actual  magnetic  curves  by  this  method  often  leads  to 
negative  initial  slopes. 

The  Rayleigh  loops,  based  on  fitting  the  end  points  and  the  initial  slope,  may  yield  good 
approximations  for  low  values  of  maximum  flux  density.  They  have  no  inflection  point, 
and  are  therefore  inherently  incapable  of  representing  a  hysteresis  curve  for  which  the 
maximum  flux  density  even  begins  to  approach  the  knee  of  the  saturation  curve.  When 
it  became  apparent  that  flux  densities  beyond  the  inflection  point  would  be  encountered, 
the  parabolas  had  to  be  abandoned  because  of  this  inherent  limitation. 
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Figure  4-19. 


4. 3. 3. 2  Cubic  Equations 

Only  cubics  wherein  B  is  a  third  degree  polynomial  in  H  can  be  considered.  The  use 
of  functions  wherein  H  is  a  third  degree  function  of  B  might  appear  attractive  at  first 
glance,  but  in  some  cases  B  would  be  a  multiple -valued  function  of  H.  Such  a  function 
would  be  an  unacceptable  representation  of  the  physical  curve. 


I 

j 


i 

I 


When  B  is  a  third  degree  function  of  H,  there  are  four  coefficients  to  be  chosen.  These 
may  be  fitted  to  the  two  end  points,  the  initial  slope,  and  the  height,  w,  of  the  loop  formed 
from  two  symmetrical  curves.  This  is  illustrated  in  Figure  4-19.  The  upper  curve  is 
obtained  by  rotating  the  bottom  one  180  degrees  about  the  origin.  The  equation  for  the 


curves  is 

B  =  ±  W 


i 

i 


where  k  is  the  initial  slope  of  either  curve,  and  the  upper  sign  is  used  for  the  upper 
curve,  and  vice  versa. 


The  above  curve  has  an  inflection  point  whose  location  depends  upon  the  curve  parameters 
w,  Pr,  Hm,  and  Bm.  Fitting  the  inflection  point  would  generally  require  a  fifth  coefficient, 
and  therefore  a  fourth  degree  polynomial. 


When  the  cubic  fit  to  actual  data  was  attempted,  the  derived  equation  was  not  monotonic. 
Therefore,  this  function  was  abandoned.  The  cubic  function  might  be  good  for  moderate 
flux  densities. 
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4. 3. 3. 3  Hyperbolic  Tangent  Function 

The  hyperbolic  tangent  functions  considered  were  of  the  form, 


B  =  Aj  +  A2  tanh  Ag  (H  +  A4) 


Such  functions  have  this  distinct  disadvantage  of  being  symmetrical  about  the  inflection 
point  H  =  -A^,  B  =  Aj.  Generally,  no  such  symmetry  exists  in  physical  curves.  The 
use  of  multiples  of  the  argument  or  powers  of  the  hyperbolic  tangent  function  will  make 
the  curve  unsymmetrical,  but  at  the  cost  of  considerable  complexity.  The  use  of  two 
different  sets  of  coefficients,  A^,  for  the  portions  of  the  curve  above  and  below  the 
inflection  point  will  also  make  the  curve  unsymmetrical.  However,  it  was  felt  that  piece- 
wise  continuous  functions  are  undesirable  and  to  be  used  only  as  a  last  resort.  If  they 
are  used,  simpler  functions  than  the  hyperbolic  tangent  should  prove  to  be  satisfactory. 
All  of  these  considerations  led  to  the  abandonment  of  the  hyperbolic  tangent  function. 


4. 3. 3. 4  Inverse  Tangent  Function 

The  inverse  tangent  function  attempted  was  of  the  form, 

tan-1  (p 

B  V 


IKstafjf-) 
m  ' 


m 


tanrK 


This  yielded  a  poor  fit  to  an  actual  curve,  in  the  vicinity  of  the  inflection  point.  The 
fit  to  actual  data  required  a  trial  and  error  process.  Furthermore,  the  function  is 
symmetrical  about  its  inflection  point,  thus  suffering  from  all  of  the  shortcomings  of  the 
hyperbolic  tangent  function.  For  these  reasons,  this  function  was  abandoned. 

4. 3. 3. 5  Froelich's  Equation 

Froelich's  equation  was  presented  in  Reference  3,  and  is  discussed  briefly  in  Reference 
4.  The  equation  is 


The  coefficients  a  and  b  were  fitted  to  one  of  the  Allegheny  Ludlum  curves.  The  two 
end  points  were  used  to  determine  the  coefficients.  Then  the  midway  point  was  tested, 
but  the  error  was  more  than  ten  percent.  Therefore,  this  form  was  abandoned.  The 
form  is  inherently  limited  to  either  concave  downward  or  concave  upward  curves.  There¬ 
fore,  at  best,  the  use  of  this  form  would  result  in  a  piecewise  continuous  function. 
Furthermore,  care  would  have  to  be  exercised  that  the  denominator  did  not  become 
infinite  within  the  H  range  used. 

4. 3. 3. 6  First  Modification  of  Froelich’s  Equation 

The  expression  considered  is 
a(H±HJ 
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where  Hc  is  the  coercive  force.  The  plus  sign  is  used  for  the  upper  curve,  and  the 
minus  sign  for  the  lower  curve. 

If  a  and  b  are  chosen  to  make  the  upper  curve  pass  through  the  points  corresponding  to 

the  residual  flux  density,  (0,  Br),  the  coercive  force,  (-Hc,  0),  and  the  positive  maximum, 

(H  ,  B  ),  then 
m’  m  ’ 


Unfortunately,  when  H  =  -Hm  is  substituted  into  this  expression,  the  result  is  B  ^  +Bm, 
for  Hc  considerably  less  than  Hm.  Evidently,  this  form  is  also  limited  as  to  the  range 
of  H,  and  would  necessitate  the  use  of  piecewise  continuous  functions.  Therefore,  it 
was  abandoned. 

4. 3. 3. 7  Second  Modification  of  Froelich’s  Equation 

The  near  successes  obtained  with  forms  similar  to  Froelich's  equation  led  to  further 
attempts  with  similar  expressions.  The  following  expression  uses  a  quadratic  denom¬ 
inator  to  better  fit  the  curve  shape: 

C3 

B  =  C4  +  - ^ . 

*  i  +  cx  h  +  c2  ir 

This  expression  was  fitted  to  most  of  the  Allegheny  Ludlum  curves  by  means  of  an 
IBM  7094  digital  computer  program.  The  program  solved  for  the  coefficients  from  four 
input  BH  points,  recalculated  these  points  to  determine  that  the  errors  in  fitting  were 
only  small  errors  due  to  roundoff,  calculated  B  for  three  interpolated  and  two  extrapolated 
values  of  H,  and  printed  out  the  real  zeroes  of  the  denominator  or  indicated  that  these 
were  complex  numbers.  When  the  zeroes  are  real,  there  is  an  infinite  discontinuity 
in  the  BH  curve.  When  the  zeroes  are  complex,  there  is  a  finite  peak  in  the  curve.  If 
either  the  peak  or  the  discontinuity  occurs  within  the  range  of  H  for  the  curve  being 
fitted,  the  expression  is  useless.  Unfortunately,  this  was  the  case  for  a  considerable 
number  of  the  Allegheny  Ludlum  curves.  Therefore,  this  functional  form  had  to  be 
abandoned,  in  spite  of  the  considerable  effort  that  had  been  expended  on  it. 
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APPENDIX  I. 

EDDY-CURRENT  LOSSES  IN  RODS 

It  will  be  shown  that  the  torque  due  to  eddy  currents  in  the  magnetic  rods  is  negligible  in 
comparison  with  that  due  to  hysteresis.  Since  the  effect  is  to  be  shown  negligible,  a 
simplified  worst-case  analysis  is  sufficient. 

Assume  that  one  rod  is  rotating  at  angular  velocity  ^  in  a  uniform  magnetic  field,  Hm, 
in  such  a  way  that  the  rod  lines  up  witn  the  lield  at  two  instants  during  one  revolution. 
The  flux  density  when  so  aligned  is  Bm.  Assuming  a  linear  relation  between  H  and  B 
as  a  gross  approximation,  the  flux  density  in  the  rod  is 

B  =  B„,  sin  o-'t.  (1-1) 

m  ' 


Its  time  derivative  is 


B  =  u)  B 

m 


cos  wt. 


(1-2) 


The  changing  flux  in  the  rod  induces  eddy  currents.  These,  in  turn,  cause  a  torque 
which  opposes  the  motion  which  caused  the  change  in  the  first  place,  in  accordance 
with  Lena's  Law.  It  is  easily  shown  that  the  torque  is 


T 


v  R4  L  B2 
8  P  u: 


(1-3) 


where  R  is  the  radius,  L  is  the  length,  and  p  is  the  electrical  resistivity  of  the  rod. 
From  Equations  (I- 1)  and  (1-2),  the  peak  torque  is 


ff  R4  L  w  B  2 
m 


VjT 


(1-4) 


For  the  rods  used,  R  is  0. 14  cm. ,  L  is  152.  5  cm. ,  and  p  is  50,  000  abohm-cm.  It 
will  be  assumed  that  u:  is  twice  the  orbital  angular  velocity,  or  0. 0C203  radians  per 
second.  The  peak  flux  density,  Bm,  is  9,  200  gauss.  The  peak  torque  is  then  0. 080 
dyne -centimeters  per  rod. 


The  peak  value  of  the  hysteresis  torque  is  calculated  for  comparison.  The  magnetic 
moment  of  one  rod  is 


M  = 


0.73 
4  IT 


Bc  = 


0. 73  R*  L 


B. 


(1-5) 


where  V  is  the  volume  of  the  rod,  and  Bc  is  the  flux  density  at  the  center.  The  factor 
0. 73  accounts  for  the  effect  of  the  gradual  reduction  of  flux  from  the  center  to  each  end 
of  the  rod.  For  the  rods  used,  which  have  a  peak  flux  density  of  9, 200  gauss,  the  peak 
value  of  the  magnetic  moment  of  one  rod  is  5,020  dyne -centimeters  per  oersted.  This 
occurs  for  an  applied  field  of  0.592  oersted.  The  peak  torque  is  therefore  2,980  dyne- 
centimeters.  This  is  more  than  four  orders  of  magnitude  greater  than  the  peak  value 
of  the  eddy-current  torque.  Therefore,  the  latter  is  negligible. 
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APPENDIX  H. 

PROGRAM  FOR  FITTING  SEVENTH-DEGREE  POLYNOMIAL 


An  IBM  7094  digital  computer  program  was  used  to  derive  a  least -squares,  seventh- 
degree  polynomial  for  each  of  the  reference  BH  curves.  The  program  will  (1)  accept 
raw  data,  that  is,  values  of  flux  density  and  internal  field  of  a  ring  sample;  (2)  convert 
to  the  external  field  applied  to  a  rod;  (3)  shift  the  data  to  an  arbitrary  origin;  (4)  solve 
the  least-squares  equations;  (5)  substitute  each  abscissa  back  into  the  derived  polynomial 
and  evaluate  it;  (6)  compute  the  error  in  the  fit  at  each  such  point;  and  (7)  printout  all 
of  these  results.  The  polynomial  is  forced  to  go  through  the  point  (0,0)  by  omitting  the 
constant  term  from  the  general  expression. 

The  input  data  are  values  of  H  and  B.  Assume  that  N  pairs  are  given: 


H, 


H. 


H 


N 


(II-l) 


Two  factors  are  also  used.  The  first  is  the  ascending-descending  index,  D^,  which  is 
+1  for  an  ascending  curve  and  -1  for  a  descending  curve.  The  polynomial  fit  is  for  an 
ascending  curve  starting  at  (0, 0).  If  the  curve  is  descending,  the  algebraic  signs  of  all 
of  the  data  values  are  reversed. 


The  second  factor  is  the  demagnetization  factor,  DM>  This  is  used  to  determine  the 
external  field  of  a  rod  corresponding  to  the  internal  field  of  a  ring  sample.  It  is  assumed 
that  the  latter  is  input,  and  the  former  is  then  computed.  In  case  the  input  is  the  ex¬ 
ternal  field,  the  factor  Dj^  is  simply  set  equal  to  zero. 

The  external  field  value  is  computed  for  each  data  point, 


H 

H 

H 


Ell^M. 
.Ej_:-Hj.:_DM 
EN  =  HN  +  °M 


(H-2) 


The  increments  for  translating  the  first  point  to  (0,  0)  are  computed.  For  the  field, 


dh  -  -  da  het 


(H-3) 


For  the  flux  density, 

db  =  -  da  bi- 

These  increments  are  applied  to  all  of  the  input  data.  For  the  field, 


(H-4) 


HS2  -  DII  +  DA  HE2» 


t 


HSj  ■  °H  +  DA  HEj’ 

— - . — 

hsn  =  dh  +  da  hen*  *n‘5) 

For  the  flux  density, 

BS2  =  DB  +  °A  B2» 

. 

BSj  =  DB  +  °A  Bj’ 

. . 

bsn  =  db  +  da  bn-  (ii_6)  I 

The  (Hg 2 }  Bgg) . (HgN  BSN^  cons^tute  N-l  pairs  of  data  to  be  used  to  fit  the 

polynomial, 

Bg  =  AjHg  A2Hg  + . +  A7HS7»  (II-7) 

I 

by  the  method  of  least  squares.  The  equations  for  the  coefficients  are 

A1  f  »srA2  f  HSf*A3  f  HsrA4  f  "sj  ‘  A5  f  "s*  +  A6  E  HSj7*A7  E  «SJ  “ f  BSjHSj- 

J  J  J  J  J  J  JJ 

i. 

A1  =  HsrA2 E  «srAs  E  HSj*A4  1  «sr  a5  e  hs’+a6  S  HS*  +  A7  E  HSi  =  E  BSiHSj2' 

J  J  J  J  J  JJJ 

I 

A1  f  Hsr  A2  f  HSj5  +  A3  f  Hsr  A4  f  HsJ  +A5  f  «S)8  +  A8  f  «S®*A7  ?  Hsl°=  E  BSjHS2>  j 

j! 

A1  f«Sj  tA2  f  »sr  A3  f«s’  +  A4  f  Hsr  A5  BHsf  *A6  BHSj10  +  A7  =  HSj  '  =  f  BSjHSj‘> 

ft 

A1  EHsrA2  SHs’*A3  Shs®  + A4  LHS2  +  A5  SHS)I0  +  A6  £HSjU+A7  EHS,12  -  iBSjHSj> 

j  j  j  j  j  j  J  J 

Al?HsrA2?HSj*A3?HS®tA4?Hs}0*A5?HSjn*A6?HS]2  +  A7?HSjS  =  ?BSJHS® 
j  j  j  j  j  j  jj 

A1  f  "S8  *  A2  E  HS)  *  A3  E  HSj10  *  A4  E  Hs"  *  A5  E  Hs}2  + A6  f  "sf  +  A7  \  HSj ' ‘  =  S  BSj  «Sj  ■ 

m-8) 

where  j  runs  from  2  to  N. 

These  are  solved  for  the  coefficients,  A.,  by  a  subroutine  used  for  simultaneous,  linear 
equations. 


11-2 


Each  abscissa,  Hgj,  is  substituted  in  turn  into  Equation  (II-7)  and  the  corresponding  value 
of  flux  density,  BCj,  is  computed, 

2  7 

BC2  =  A1  ^S2  +  A2  HS2  +  •  •  •  •  +  A7  HS2  ’ 

BCj  =  A1  HSj  +  A2  HSj2  + . +  A7  HSj7» 

2  7 

BCN  =  A1 HSN  +  A2  HSN  +  •  •  •  •  +  A7  hsn  ‘  (n_9) 

The  error  in  the  polynomial  at  each  point  is  computed.  It  is  the  difference  between  the 
flux  density  computed  from  the  polynomial  and  the  corresponding  input  value, 

BE2  =  BC2  ‘  BS2’ 

BEj  =  BCj  ‘  BSj» 

BEN  =  BCN  ‘  BSN'  (H-10) 

The  computed  values  of  the  flux  density  are  retranslated,  by  subtracting  DB> 

BR1  =  ■  DB’ 

BR2  =  BC2  “  °B’ 

BRj  =  BCj  "  DB’ 

brn=bcn~db*  (n-ii) 

The  argument  of  the  fitted  polynomial  remains  the  translated  value  of  the  field,  which  is 
zero  at  the  beginning  point  of  the  curve.  This  is  the  same  way  that  the  polynomial  is 
used  in  the  main  program. 

It  should  be  emphasized  that  the  fitted  polynomial  pertains  to  an  ascending  curve. 
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APPENDIX  HI. 

ELEVEN-POINT  FIT  OF  THE  RATIONAL  FRACTION 


This  appendix  lists  the  equations,  input,  and  output  of  an  IBM  7094  program  used  for 
deriving  the  coefficients,  A^,  of  the  function 


B  =  L 


E  Ai 
i  (H-Hj)2  +  C2 


(m-i) 


This  function  may  be  used  as  an  analytical  approximation  for  a  BH  curve,  as  discussed 
in  Section  4. 3. 2. 1. 

The  range  of  H  is  divided  by  20,  to  obtain  21  points.  The  corresponding  21  values  of  B, 
from  the  curve  to  be  fitted,  are  also  tabulated.  The  eleven  alternate  pairs  of  values, 
beginning  with  the  first  pair,  are  used  for  fitting.  The  intermediate  ten  pairs  of  values 
are  used  for  checking  the  goodness  of  the  fit. 

A  rearrangement  of  Equation  (III-l)  is  convenient.  First,  the  increment  between  con¬ 
secutive  values  of  H  (1/20  of  the  range)  is  designated  AH.  Also,  C  is  expressed  as  a 
multiplier,  H^,  times  AH.  Then 
o  A. 


B(AH)2  =  E  -s-L 


i  n  +  H. 


(ID- 2) 


where  n  is  zero  or  a  positive  integer.  The  method  of  summing  is  explained  later. 

The  program  may  be  used  for  fitting  either  ascending  or  descending  magnetic  hysteresis 
curves.  An  index,  D^,  i3  set  equal  to  +1  if  the  curve  is  ascending  and  -1  if  it  is  descending. 

The  input  to  the  program  consists  of  the  index,  the  multiplier,  H^,  and  the  21  pairs 
of  values : 

Ho'  Bo. 

Hj,  Bj, 


H20»  B20* 


(m-3) 


The  first  computations  are  all  of  the  denominator  factors  which  will  be  required.  These 
are  for 

n  =  0,  1,  2, . 20.  (m-4) 


III-l 


The  factors  are 


H0  =  0  * 

0  H* 


R,  = 


i2-„k2  ’ 


Ro  = 


1 


2  9  > 

2  +  hk 


_ 1 

on2  u  2  * 

20  +  Hj£ 


I 


(m-5) 


The  coordinates  of  the  eleven  points  to  be  fitted  are  substituted  into  Equation  (III-2). 

For  any  general  point,  one  of  the  terms  will  correspond  to  n  =  0,  the  two  adjacent  terms 
will  correspond  to  n  =  1,  etc.  The  scheme  is  apparent  from  an  examination  of  the  sub¬ 
scripts  in  the  determinant  of  the  system  of  equations.  Since  alternate  points  are  used, 
only  R  factors  with  even  subscripts  appear.  The  eleven  simultaneous  linear  equations 
for  the  Aj  are,  in  matrix  form: 


Ro 

R2 

R4 

R6 

R8 

R10 

R12 

R14 

R16 

Rl8 

R2 

Ro 

R2 

R4 

R6 

R8 

R10 

R12 

R14 

R16 

R4 

R2 

Ro 

R2 

R4 

R6 

R8 

R10 

R12 

R14 

R6 

R4 

R2 

Ro 

R2 

R4 

R6 

R8 

R10 

R12 

R8 

R6 

R4 

R2 

Ro 

R2 

R4 

R6 

R8 

R10 

R10 

R8 

R6 

R4 

R2 

Ro 

R2 

R4 

R6 

R8 

R12 

R10 

R8 

R6 

R4 

R2 

Ro 

R2 

R4 

R6 

R14 

R12 

R10 

R8 

R6 

R4 

R2 

Ro 

R2 

R4 

R16 

R14 

R12 

R10 

R8 

R6 

R4 

R2 

Ro 

R2 

R18 

R16 

R14 

R12 

R10 

R8 

R6 

R4 

R2 

Ro 

R20 

R18 

R16 

R14 

R12 

R10 

R8 

R6 

R4 

R2 

R20 

Ao 

Bo 

R18 

A1 

B2 

R16 

A2 

B4 

R14 

A3 

B6 

R12 

A4 

B8 

R10 

A5 

=  Da(AH)2 

B10 

R8 

A6 

B12 

R6 

A7 

B14 

R4 

A8 

B16 

R2 

A9 

B18 

Ro 

A10 

B20 

(HI-6) 


After  solving  for  the  eleven  A^,  the  computed  and  the  intermediate  points  are  tested, 
using  Equation  (HI-2)  solved  for  B.  The  values  of  B  so  computed  are  designated  B^. 
The  twenty-one  B^'s  are  computed  from  the  following  equations: 


HI-2 


(AH)' 


(AH) 


(AH)' 


T  <RoAo  +  R2A1  +  R4A2  +  R6A3  +  R8A4  +  RlOA5  +  R12A6 
+  R14A7  +  R16Ag  +  RjgAg  +  R20Ai0^ 

2  (RlAo  +  R1A1  +  R3A2  +  R5A3  +  R7A4  +  R9A5  +  R11A6 
+  Rl3A7  +  Rl5Ag  +  Rl7Ag  +  R^qAjq), 

1  (Vo  +  RoAl  +  R2A2  *  R4A3  +  V4  +  R8A5 
+  **1(^6  +  R12A7  +  R14A8  +  R16A9  +  Rl8AlC^> 


(AH) 


2  (R3Ao  +  R1A1  +  R1A2  +  R3A3  +  R5A4  +  R7A5 


+  RgAg  +  R11A?  +  Rl3A8  +  Rl5Ag  + 


(AH) 


2  (R4Ao  +  R2A1  +  RoA2  +  R2A3  +  R4A4  +  R6A5 


(AH)' 


(AH)' 


+  RgAg  +  R10A?  +  Rl2Ag  +  Rl4Ag  +  Ri6A10^ 

(R5Ao  +  R3A1  +  RlA2  +  R1A3  +  R3A4  +  R5A5 
+  R?Ag  +  RgA?  +  RuAg  +  Rl3Ag  +  Rl5Al0), 

(R6Ao  +  R4A1  +  R2A2  +  RoA3  +  R2A4  +  R4A5 

+  RgAg  +  RgA?  +  RjgAg  +  Rl2Ag  +  ri4A!0^ 


(AH) 


2  (R7Ao  +  Vl  +  R3A2  +  R1A3  +  R1A4  +  R3A5 


+  R5A6  +  V*  +  RQAR  +  RnAQ  +  RisAm)> 


9^8  +  119  +  13  10' 


(AH) 


2  (R8Ao  +  R6A1  +  R4A2  +  R2A3  +  RoA4  +  R2A5 
+  R4A6  +  R6A7  +  R8A8  +  R10A9  +  R12AlO)> 


(AH) 


2  (R9Ao  +  Vl  +  R5A2  +  R3A3  +  R1A4  +  R1A5 
+  R3A6  +  R5A7  +  W7A8  +  V9  +  R11A10)» 


m-3 


BC10  = 


BC11  = 


BC12  " 


BCl3  = 


BC14  = 


BC15  = 


BC16  = 


BC17  = 


BC18  " 


BC 19  = 


(AH)' 


D, 


(AH)' 


(AH  )' 


(AH)' 


(AH)' 


D. 


(AH)' 


<RlO*o  +  R8A1  +  R6A2  +  R4A3  +  R2A4  +  RoA5 
+  RgAg  +  R4A7  +  RgAg  +  RgAg  +  R^jA^), 

(RllAo  +  +  R7A2  +  R5A3  +  R3A4  +  R1A5 

+  R^Ag  +  RgA7  +  RgAg  +  R?Ag  +  RgAlC), 

^R12Ao  +  R10A1  +  R8A2  +  R6A3  +  R4A4  +  R2A5 
+  RoA6  +  R2A7  +  R4A8  +  R6A9  +  R8A1(P> 

^Rl3Ao  +  R11A1  +  R9A2  +  R7A3  +  R5A4  +  R3A5 

+  R^Ag  4-  R^Arj  4“  RgAg  4*  RgAg  4-  Jq)) 

^R14Ao  +  R12A1  +  R10A2  +  R8A3  +  R6A4  +  R4A5 
+  RgAg  +  R0A?  +  RgAg  +  R4Ag  +  RgAl0), 

^Rl5Ao  +  Rl3Al  +  R11A2  +  R9A3  +  R7A4  +  R5A5 

+  RgAg  +  RjA  7  +  RjAg  +  RgAg  +  R5Al0), 


^R16Ao  +  R14A1  +  R12A2  +  R10A3  +  R8A4  +  R6A5 

+  R4Ag  +  RgApj  +  RQAg  +  RgAg  +  R4AlQ), 


°A 

(ar)2  (R17Ao  +  R15A1  +  Rl3A2  +  R11A3  + 


R9A4  +  R7A5 


+  R5A6  +  R3A7  +  R1A8  +  R1A9  +  R3A10^ 


^Rl8Ao  +  R16A1  +  R14A2  +  Rl2A3  +  **1(^4  +  **8^5 
+  RgAg  +  R4A7  +  RgAg  +  RQAg  +  RgAjg), 

^R19Ao  +  R17A1  +  R15A2  +  Rl3A3  +  R11A4  +  R9A5 

+  R?Ag  4  RgA7  4-  RgAg  4-  RjAg  4-  RjAjq), 


III-4 


(m-7) 


da 

Bc20=  (Ah)2  ^R2oA°  +  RlsAl  +  R16A2  +  R14A3  +  R12^4  +  R10A5 

+  R8A6  +  R6A7  +  R4A8  +  ^2^9  +  RoAlO^ 

The  error  in  each  computed  value  is  the  difference  between  that  value  and  the  actual 
(input)  value  of  B.  These  errors  are  computed  from 

BE0  =  BC0  "  '  » 

EE1  =  BC1  ‘  Bl> 

BE2  =  BC2  ‘  B2’ 

BE20=  BC20  ‘  B20* 

The  output  includes  the  eleven  coefficients,  AQ,  A^,  A2, . A^q,  and  the  tabulated 

input  pairs  of  values,  computed  flux  density,  and  error  in  flux  density: 


Ho 

Bo 

BCo 

BEo 

H1 

B1 

BC1 

bei 

H2 

B2 

BC2 

BE2 

H20 

B20 

BC20 

BE20- 
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